信頼区間について

2023年2月26日 15:29

以前の記事で、「専門家でも分かるように信頼区間の説明をしてください」という面接質問に答えるのは難しいのではないか、という話をしました。それでちょっと昔の教科書などを引っ張り出してきて、復習がてら自分なりに整理してみました。

そういえば医者になった知人に久しぶりに会って、数学・統計やってたよ、と話をしたら「p値って何なの？」と目を輝かせて質問されました。医者でも理解は難しい概念です。

まずは例

昔の教科書を引っ張り出してきたのですが、信頼区間に関する説明があったよなー、と記憶していた講義ノートを読んでも全然出てこない。困ったなー、と思ったら以前買った本に良い例がありましたので、それを参考に例を作ってみました。

Teaching StatisticsというAndrew Gelmanの本です。教科書というよりは、どうやって統計学を大学などで教えるのか、という目的で書かれた本のようです。私が持っているのは1st editionでしたが、第二版では９章の9.4 "Confidence Intervals: theory"という項にわかりやすい説明があったので、これを踏まえて紹介します。

設定

講義に参加している学生に紙を配ります。紙には自分の体重を書いてもらいます。書いたら箱に入れて（例では帽子とありました）もらい、回収します。その間に一人の学生にボランティアとして参加してもらい、箱の中にある体重を全てメモし、平均と標準偏差を計算します。一方ほかの学生には、一人一人はこの中の紙をよく混ぜて、5枚だけ取り出し記録された体重を書き写してもらいます。５つの数字を使って平均、標準偏差を計算してもらいます。最後に、ボランティア学生が全部の紙を使って計算が終わったら、分散だけを発表してもらいます。

さらに、発表された分散を使って、各学生には68%信頼区間を計算してもらいます。

シミュレーションしてみる

実際に紙に書いてもらったり何人もの学生を集めたりはできないので、シミュレーションしてみます。また、体重は嫌がる学生もいるかもしれないので、身長という設定にしてみます。シミュレーションにはRを使います。

e-Statの「国民健康・栄養調査」に「身長・体重の平均値及び標準偏差 - 年齢階級，身長・体重別，人数，平均値，標準偏差 - 男性・女性，１歳以上〔体重は妊婦除外〕」という統計があります。これによると、2019年のデータで20歳の男性170.2、女性158.6cmとなっています。クラスにいる男性、女性が半々だと想定すると、平均は男性・女性の平均で164.4cmと考えられます。40人のクラス全員の身長のデータは、乱数を発生させて作ります。標準偏差はわからないので、大体±20cmぐらいに収まるだろうという予想の元、10に設定します。

set.seed(893)
y <- rnorm(40, 164.4, 10)

グラフを見てみると、140cm台の人が3人、180センチを超える人が一人なので、まあまあよさそうなデータが作れています。

hist(y)

全体のデータの集計

一人の生徒にボランティアで、40人分のデータを集計してもらいます。平均と標準偏差を計算してもらいます。これが母集団の、実際に推定したい平均値と標準偏差になります。平均の数字はまだ発表してもらわないで、標準偏差だけを発表してもらういます。

mean(y); sd(y)

標準偏差は10.5となりました。

残りの生徒それぞれに信頼区間を計算してもらう

一人の生徒が40人分のデータを集計する間に、残りの生徒に箱を回して、5枚紙を出して身長データをノートに記録してもらいます。5つのデータを取るには、Rだとsample関数でできます。これを39人の生徒分繰り返したいので、replicate関数を使います。

sample(y, 5, replace=FALSE) # 生徒一人分のデータ集計作業をする。５つのサンプルを取る。

# この作業を39回繰り返す。結果は行列に入る。
height_sampling <- replicate(39, sample(y, 5, replace = FALSE))

それぞれの５サンプルのデータから、平均と標準偏差を計算します。

height_mean <- apply(height_sampling, 2, mean)
height_sd <- apply(height_sampling, 2, sd)

この頃にはボランティアの生徒が40人分のデータ集計を終えていると思われるので、標準偏差を発表してもらいます。すでに計算した通り10.5が答えです。発表された標準偏差をもとに68%信頼区間を計算します。各生徒のデータは５サンプルなので、10.5を５の平方根で割ります。

height_conf_interval <- tibble(sample_id = seq(1, dim(height_sampling)[2]),
                               height_mean = height_mean,

                               height_lower_bound = height_mean - 10.5/sqrt(5),
                               height_upper_bound = height_mean + 10.5/sqrt(5))
height_conf_interval |> head()


sample_id height_mean height_lower_bound height_upper_bound height_lower_bou…¹ heigh…²
<int>       <dbl>              <dbl>              <dbl>              <dbl>   <dbl>
1         1        159.               154.               163.               156.    162.
2         2        155.               150.               160.               149.    161.
3         3        163.               158.               167.               158.    167.
4         4        164.               160.               169.               162.    167.
5         5        160.               155.               165.               155.    165.
6         6        164.               159.               169.               159.    168.

これらの信頼区間をグラフにします。グラフ化することで、それぞれの信頼区間は５つのデータを取った時に、どの５データが取れるかによってランダムに決まるということがわかります。ちなみに縦の線が、40人全員から計算してもらった真の平均値です。今回のデータだと162.4cmです。

height_conf_interval |> ggplot(aes(x=height_mean, y = sample_id)) +
  geom_pointrange(aes(xmin=height_lower_bound, 
                      xmax=height_upper_bound)) + 
  geom_vline(xintercept = mean(y)) + 
  labs(title="Confidence Intervals calculated by each student", 
       x="Height", y = "Sample ID")

これをみると、「信頼区間は68％の確率で真の平均を含む」という言い方は間違いであるとわかります。統計解析で得られる信頼区間は、上記のような信頼区間のどれか１つです。１つの信頼区間を見ると、真の平均の162.4cmを含むか含まないかのいずれかです。

95％信頼区間などで言われる９５％という確率は、求められた信頼区間に対するものではなくて、この手順で何回も信頼区間を計算すれば、68%の確率で真の平均を含む、という意味です。今回の例だと39個信頼区間を作ってもらったので、大体27個（36 x 68%）の信頼区間が真の平均を含むということが言えます。

ちなみに真の平均162.4を含むのは

height_conf_interval |>
  mutate(is_true_mean_covered = ifelse(162.4>=height_lower_bound & 162.4<=height_upper_bound, 1, 0)) |>
  group_by() |>
  summarise(true_mean_covered_count = sum(is_true_mean_covered))

24

24回だったので、27回には少し足りないですが、まあまあ68％に近い回数にはなっています。

分散未知の場合

一人の生徒に40人分のデータから計算した標準偏差を発表してもらったのは、統計を学んだ人なら一度は聞いたことがある「分散既知」の場合を指しています。母集団（学生40人）の分散だけはわかっている、という状態を作るために発表してもらいました。昔統計を勉強した時には、どんな状況でそんなことが起きるのか、と疑問に思っていました。

仮に母集団の分散（標準偏差）もわからない場合、つまり一人のボランティアに標準偏差を発表してもらわない場合は、各学生は５つのサンプルを撮ったので、自由度４のt分布を使って信頼区間を計算します。

height_conf_interval <- height_conf_interval |>
    mutate(height_lower_bound_t = height_mean - qt(0.84, 4, 0)*height_sd/sqrt(5),
         height_upper_bound_t = height_mean + qt(0.84, 4, 0)*height_sd/sqrt(5))

パッと見たところそんなに変わらないですが、少しだけ信頼区間の幅が広くなっています。

まとめ

ということで信頼区間の復習がてら、"Teaching Statistics"の例を使いながら紹介しました。