R言語で分散が既知のときの母平均の信頼区間を算出する関数やその実行例について解説します。
実行例では、母平均の信頼区間だけでなく母平均の差の信頼区間の求め方について見ていきます。母平均の差の信頼区間の実行れでは、対応のあるデータの母平均の差の信頼区間の計算の仕方も紹介します。
この記事で扱うプログラミングコードは以下からダウンロードできます。
母平均および母平均の差の信頼区間差については、次の記事を参照してください。
-
【統計学】母平均の信頼区間
標本が正規分布に従う場合の母平均の信頼区間の導出についてみていきます。 母分散が既知または未知である場合についてそれぞれの信頼区間を解説します。 母平均の差の信頼区間については以下を参照。 信頼区間 ...
続きを見る
-
【統計学】母平均の差の信頼区間
標本が正規分布に従う場合の母平均の差の信頼区間の導出についてみていきます。 母分散が既知または未知である場合について二群の母平均の差の信頼区間を解説します。 母平均の信頼区間については以下の記事を参照 ...
続きを見る
分散が既知のときの母平均の信頼区間
分散が既知のときの母平均の信頼区間とR言語の母平均の信頼区間の関数についてまとめました。
母平均と母平均の差の信頼区間の概要について列挙しています。
母平均の信頼区間
\(x_1, x_2, \ldots, x_n\)を\(N(\mu, \sigma^2)\)からの\(n\)個の無作為標本とする。また、標本平均と不偏標本分散を次とする。
母平均の信頼区間
母平均\(\mu\)の\(100(1 - \alpha)\)%信頼区間は次で与えられる。
母平均の差の信頼区間
\(x_{11}, x_{12}, \ldots, x_{1n_1}\)を\(N(\mu_1, \sigma_1^2)\)からの\(n_1\)個の無作為標本、\(x_{21}, x_{22}, \ldots, x_{2n_2}\)を\(N(\mu_2, \sigma_2^2)\)からの\(n_2\)個の無作為標本とする。
母平均の差の信頼区間
対応のないデータかつ2母集団が独立のとき
2つの)のとき、母平均の差\(\mu_1 - \mu_2\)の\(100( 1 - \alpha)\)%信頼区間は次で与えられる。
対応のあるデータのとき
2変量正規母集団\(N(\boldsymbol{\mu}, \boldsymbol{\Sigma})\)からの大きさ\(n\)の無作為標本\((x_{11}, x_{21}), \ldots, (x_{1n}, x_{2n})\)が得られたとする。ここに
共分散行列\(\boldsymbol{\Sigma}\)が既知であり対応のあるとき、母平均の差\(\mu_1 - \mu_2\)の\(100(1 - \alpha)\)%信頼区間は次で与えられる。
ここに\(n_1 = n_2 = n\)、\(x_{di} = x_{1i} - x_{2i},\ i = 1, \ldots, n\)、
また、\(\sigma_d^2\)は\(X_{di} = X_{1i} - X_{2i},\ i = 1,\ldots, n\)の分散\(\sigma_d^2 = \mathrm{Var}[X_{di} ]= \sigma_1^2 + \sigma_2^2 - 2\sigma_{12}^2\)である。
関数z.test
パッケージBSDAの関数z.testで母平均および母平均の差の信頼区間を計算することができます。以下、z.testとその引数です。
| x | numeric型のベクトル。標本（観測値）。NAまたはInfは除去される。 |
| y | numeric型のベクトル。別のグループの標本（観測値）。NAまたはInfは除去される。 |
| alternative | character型。"two.sided"で両側検定、"greater"で右片側検定、"less"で左片側検定を設定する。 |
| mu | numeric型。帰無仮説の下での母平均の値、または母平均の差の値。 |
| sigma.x | numeric型。xの標準偏差。 |
| sigma.y | numeric型。yの標準偏差。 |
| conf.level | numeric型。信頼区間を計算する際の信頼係数。 |
実行例
母平均の信頼区間の実行例について紹介します。
母平均の信頼区間の算出のためにパッケージBSDA、グラフの作成の際にパッケージggplot2を用います。あらかじめインストールしておいてください。
1 2 3 4 | install.packages("ggplot2") install.packages("BSDA") library(ggplot2) library(BSDA) |
母平均の信頼区間
次の木に関するデータセットtreesを例に母平均の信頼区間の計算を行っていきます。
1 2 3 4 5 6 7 8 9 | > dataset <- trees > head(dataset) Girth Height Volume 1 8.3 70 10.3 2 8.6 65 10.3 3 8.8 63 10.2 4 10.5 72 16.4 5 10.7 81 18.8 6 10.8 83 19.7 |
Heightの95%信頼区間を算出するには、次のように関数z.testの引数conf.levelを0.95、sigma.xに分散を与えて後ろに$conf.intを付けます。
分散を40にしているのは、標本分散がだいたい40ぐらいだからです。
1 2 3 | sigma_squared <- 40 conf_level <- 0.95 ci <- z.test(dataset$Height, sigma.x = sqrt(sigma_squared), conf.level = conf_level)$conf.int |
ciを参照すると信頼区間の下限と上限、また信頼水準が格納されているのが分かります。
1 2 3 4 | > ci [1] -7.5373939 0.1373939 attr(,"conf.level") [1] 0.95 |
信頼区間の上限と下限を取得するには、次のようにciの1番目と2番目をインデックスで参照すれば大丈夫です。
1 2 | ci_lower <- ci[1] #信頼区間の下限 ci_upper <- ci[2] #信頼区間の上限 |
計算した信頼区間を他のファイルに保存したいときは、次のように関数data.frameでデータフレームに計算した値を格納しwrite.csvなどで他のファイルへ保存します。
1 2 3 | resultTable <- data.frame(sample.mean = mean(dataset$Height), 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 | ggplot(dataset, aes(x = 1, y = Height)) + geom_dotplot(binaxis = "y", stackdir = "center") + stat_summary(aes(x = 1, y = Height), fun = mean, fun.min = function(x) mean(x), fun.max = function(x) mean(x), geom = "crossbar", width = 0.3, color = "#F8766D") + stat_summary(aes(x = 1, y = Height), fun.min = function(x) ci[1], fun.max = function(x) ci[2], geom = 'errorbar', width = 0.25, color = "#00BFC4", size = 1) + geom_text(data = data.frame(x = c(mean(dataset$Height), ci), name = factor(c("sample mean", paste0(100 * conf_level, rep("% c.i.", 2))), levels = c("sample mean", paste0(100 * conf_level, "% c.i.")))), aes(x = 0.7, y = x, color = name, label = paste0(c("sample mean", "c.i.lower", "c.i.upper"), ": ", round(x, 3)))) + labs(color = "") |
母平均の信頼区間
- 赤線はHeightの標本平均
- 青線はHeightの母平均の上限と下限
を表しています。
母平均はだいたい73から78の範囲にあることを推定できました。
母平均は母集団のパラメータであり、上の図の赤線の標本平均ではないことに注意してください。
母平均の差の信頼区間（2つの母集団が独立で対応のない場合）
続いて母平均の差の信頼区間の算出例について見ていきます。
次のデータセットToothGrowthのsuppがVCであるグループとOJであるグループのLenの母平均の差の信頼区間を求めたいと思います。
1 2 3 4 5 6 7 8 9 | > dataset <- ToothGrowth > head(dataset) len supp dose 1 4.2 VC 0.5 2 11.5 VC 0.5 3 7.3 VC 0.5 4 5.8 VC 0.5 5 6.4 VC 0.5 6 10.0 VC 0.5 |
母平均の差の信頼区間を計算するには、関数z.testの引数を次のように与えます。
1 2 3 4 5 6 7 8 9 | sigma_x_squared <- 70 sigma_y_squared <- 45 conf_level <- 0.95 len_VC <- dataset$len[dataset$supp == "VC"] len_OJ <- dataset$len[dataset$supp == "OJ"] ci <- z.test(len_VC, len_OJ, sigma.x = sqrt(sigma_x_squared), sigma.y = sqrt(sigma_y_squared), conf.level = conf_level)$conf.int |
suppがVCであるグループとOJであるグループの母平均の差の信頼区間は次の画像の通りです。
1 2 3 | resultTable <- data.frame(diff = mean(len_VC) - mean(len_OJ), c.i.lower = ci[1], c.i.upper = ci[2], row.names = NULL) write.csv(resultTable, "母平均の差の信頼区間_分散既知.csv", row.names = FALSE) |
母平均の差の信頼区間
2つのグループの信頼区間をプロットする際は、次のようにggplotを実行します。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 | ci_OJ <- z.test(len_OJ, sigma.x = sqrt(sigma_x_squared), conf.level = conf_level)$conf.int ci_VC <- z.test(len_VC, sigma.x = sqrt(sigma_y_squared), conf.level = conf_level)$conf.int dataset_ci <- data.frame(mean_len = tapply(dataset$len, dataset$supp, mean), supp = as.factor(c("OJ", "VC")), ci = c(dist(ci_OJ), dist(ci_VC)) / 2) ggplot() + geom_dotplot(data = dataset, aes(x = supp, y = len, color = supp, fill = supp), binaxis = "y", stackdir = "center") + geom_line(data = dataset_ci, aes(x = supp, y = mean_len, group = 1), size = 1, alpha = 0.5) + geom_point(data = dataset_ci, aes(x = supp, y = mean_len, group = 1), size = 2, alpha = 0.5) + geom_errorbar(data = dataset_ci, aes(x = supp, y = mean_len, ymin = mean_len - ci, ymax = mean_len + ci), width = 0.25, size = 1, alpha = 0.5) + geom_text(data = dataset_ci, aes(x = supp, y = mean_len - ci + 1, label = paste0(conf_level, "% c.i.lower: ", round(mean_len - ci, 3)))) + geom_text(data = dataset_ci, aes(x = supp, y = mean_len + ci + 1, label = paste0(conf_level, "% c.i.upper: ", round(mean_len + ci, 3)))) + geom_text(data = dataset_ci, aes(x = supp, y = mean_len + 1, label = paste0("Sample mean: ", round(mean_len, 3)))) |
母平均の差の信頼区間
上の画像から分かるように、赤のOJのグループの上限と青のVCのグループの下限の差が母平均の差の信頼区間の上限、OJのグループの下限と青のVCのグループの上限の差が母平均の差の信頼区間の下限に近い値をとっているのが分かります。
各グループの母平均の信頼区間を描画したものであり母平均の差の信頼区間のプロットでないことに注意してください。
母平均の差の信頼区間（対応のあるデータ）
最後に対応のあるデータの母平均の差の信頼区間の計算例を紹介します。
次のデータのsleepについて母平均の差の信頼区間を計算していきます。
1 2 3 4 5 6 7 8 9 | > dataset <- sleep > head(dataset) extra group ID 1 0.7 1 1 2 -1.6 1 2 3 -0.2 1 3 4 -1.2 1 4 5 -0.1 1 5 6 3.4 1 6 |
対応のあるデータの母平均の差の信頼区間を計算するには、関数z.testの引数を次のように与えます。
ポイント
対応のあるデータの母平均の差の信頼区間を計算したい場合は、次のように差をとったデータを引数xに渡し、標準偏差の引数sigma.xを「第一母集団の分散+第二母集団の分散-2×2つの母集団の共分散」にします。
1 2 3 4 5 6 7 8 | conf_level <- 0.95 Sigma <- matrix(c(3, 3, 3, 4), nrow = 2) extra1 <- dataset$extra[dataset$group == 1] extra2 <- dataset$extra[dataset$group == 2] ci_paired <- z.test(extra1 - extra2, sigma.x = sqrt(Sigma[1, 1] + Sigma[2, 2] - 2 * Sigma[1, 2]), conf.level = conf_level)$conf.int |
次のように、母平均の信頼区間と同様に母平均の差の信頼区間の上限と下限の取得が可能です。
1 2 | ci_lower <- ci_paired[1] #信頼区間の下限 ci_upper <- ci_paired[2] #信頼区間の上限 |
まとめ
R言語で分散が既知のときの母平均の信頼区間を計算する関数やその実行例をみていきました。
パッケージBSDAの関数z.testを用いることで分散が既知のときの母平均およぼ母平均の差の信頼区間を計算することができます。
