R言語で母分散の信頼区間を算出する関数やその実行例について解説します。
実行例では、データセットから母分散の信頼区間を計算するだけでなく、母分散の信頼区間の意味を図で説明します。
この記事で扱うプログラミングコードは以下からダウンロードできます。
母分散の信頼区間差については、次の記事を参照してください。
母分散の信頼区間
分散の信頼区間とR言語の関数についてまとめました。
母分散の信頼区間の概要について列挙しています。
母分散の信頼区間
正規母集団\(N(\mu, \sigma^2)\)からの大きさ\(n\)の無作為標本\(x_1, \ldots, x_n\)が得られたときの、母分散\(\sigma^2\)の\(100(1- \alpha)\)%信頼区間を紹介する。以降、標本平均を\(\bar{x} = (1 /n) \sum_{i = 1}^n x_i\)、不偏標本分散を\(u^2 = \sum_{i=1}^n (x_i - \bar{x})^2 / (n-1)\)、自由度\(n\)のカイ二乗分布の上側\(\alpha\)点を\(\chi_{n, \alpha}^2\)とする。
母分散の信頼区間
母分散\(\sigma^2\)の\(100(1 - \alpha)\)%信頼区間は次で与えられる。
関数varTest
パッケージEnvStatsの関数varTestを用いることで母分散の検定を行うことができます。
以下varTestとその引数です。
varTest(x, alternative = "two.sided", conf.level = 0.95, sigma.squared = 1, data.name = NULL)
| 引数 | 説明 |
| x | numeric型のベクトル（観測値）。NAやNaN、infにも適用できるが除外される。 |
| alternative | character型。"two.sided"で両側検定、"greater"で右片側検定、"less"左片側検定を指定できる。 |
| conf.level | 信頼区間の信頼度数。 |
| sigma.squared | 帰無仮説のときの母分散の値。 |
| data.name | character型。検定に用いるデータの名前。 |
実行例
母分散の信頼区間の実行例について紹介します。
母分散の信頼区間の算出のためにパッケージEnvStats、グラフの作成の際にパッケージggplot2、gridExtraを用います。あらかじめインストールしておいてください。
1 2 3 4 | install.packages("ggplot2") install.packages("EnvStats") library(ggplot2) 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 |
SpeciesがsetosaであるSepal.Lengthのヒストグラムは以下の通りです。
1 2 3 4 5 6 7 8 9 10 11 | boxPlot <- ggplot(dataset[dataset$Species == "setosa", ], aes(x = 1, y = Sepal.Length)) + geom_boxplot() + scale_y_continuous() + coord_flip() + geom_point(position = position_dodge(width = 0.9), alpha = 0.3, size = 1) histPlot <- ggplot(dataset[dataset$Species == "setosa", ], aes(x = Sepal.Length)) + geom_histogram(bins = 15) + scale_y_continuous() grid.arrange(histPlot, boxPlot) |
SpeciesがsetosaであるSepal.Lengthのヒストグラム
SpeciesがsetosaであるSepal.Lengthの母分散の信頼区間の計算を行いたいと思います。
setosaに限定する理由は、すべてのSpeciesを含めてしまうと正規性が成り立たなくなるためです。ヒストグラムを描けば分かりますが、全てを含めた場合Sepal.Lengthのヒストグラムは多峰性を持ちます。
SpeciesがsetosaであるSepal.Lengthの95%信頼区間を算出するには、次のように関数varTestの引数conf.levelを0.95を与えて後ろに$conf.intを付けます。
1 2 3 | conf_level <- 0.95 ci <- varTest(dataset$Sepal.Length[dataset$Species == "setosa"], conf.level = conf_level)$conf.int |
ciを参照すると信頼区間の下限と上限、また信頼水準が格納されているのが分かります。
1 2 3 4 5 | > ci LCL UCL 0.08669881 0.19293982 attr(,"conf.level") [1] 0.95 |
信頼区間の上限と下限を取得するには、次のようにciの要素を名前で取得すれば大丈夫です。
1 2 | ci_lower <- ci["LCL"] #信頼区間の下限 ci_upper <- ci["UCL"] #信頼区間の上限 |
計算した信頼区間を他のファイルに保存したいときは、次のように関数data.frameでデータフレームに計算した値を格納しwrite.csvなどで他のファイルへ保存します。
1 2 3 | resultTable <- data.frame(sample.variance = var(dataset$Sepal.Length[dataset$Species == "setosa"]), c.i.lower = ci["LCL"], c.i.upper = ci["UCL"], 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 | dataset_ci <- data.frame(x = c(mean(dataset$Sepal.Length[dataset$Species == "setosa"]) + c(sd(dataset$Sepal.Length[dataset$Species == "setosa"]), -sd(dataset$Sepal.Length[dataset$Species == "setosa"])), mean(dataset$Sepal.Length[dataset$Species == "setosa"]) + c(sqrt(ci), -sqrt(ci))), name = factor(c(rep("SD", 2), paste0(100 * conf_level, paste0("% c.i.", rep(c("upper", "lower"), 2)))), levels = c("SD", paste0(100 * conf_level, "% c.i.", c("upper", "lower"))))) ggplot(dataset[dataset$Species == "setosa", ], aes(x = 1, y = Sepal.Length)) + geom_dotplot(binaxis = "y", stackdir = "center") + stat_summary(aes(x = 1, 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 = 0.95, y = Sepal.Length), fun.min = function(x) mean(x) - sqrt(ci["LCL"]), fun.max = function(x) mean(x) + sqrt(ci["LCL"]), geom = 'errorbar', width = 0.1, color = "#619CFF", size = 1) + stat_summary(aes(x = 1.05, y = Sepal.Length), fun.min = function(x) mean(x) - sqrt(ci["UCL"]), fun.max = function(x) mean(x) + sqrt(ci["UCL"]), geom = 'errorbar', width = 0.1, color = "#F8766D", size = 1) + geom_text(data = dataset_ci, aes(x = 0.7, y = x, color = name, label = paste0(c(rep("SD", 2), rep(c("c.i.lower", "c.i.upper"), 2)), ": ", round(x, 3)))) + labs(color = "") + scale_colour_manual(values = c("#00BA38", "#619CFF", "#F8766D")) |
母分散の信頼区間の意味
- 緑線はsetosa種のSepal.Lengthの標本平均±標準偏差
- 赤線はsetosa種のSepal.Lengthの標本平均±標準偏差の信頼区間の上限
- 青線はsetosa種のSepal.Lengthの標本平均±標準偏差の信頼区間の下限
を表しています。
母分散の信頼区間がだいたい0.08から0.19であることから、その平方根である0.29から0.44が標準偏差の信頼区間となります。
データのばらつきの指標である標本平均±標準偏差は上の青から赤の範囲にあることが分かります。（標本平均は考慮していません。）
まとめ
R言語で母分散の信頼区間を計算する関数やその実行例をみていきました。
パッケージEnvStatsの関数varTestを用いることで母分散の信頼区間を計算することができます。
