第6回講義資料
統計的推定
スライド
セットアップ
今回使用するパッケージは{tidyverse}のみだ。データはinference_data.csvであり、LMSから入手できる。読み込んだデータはdfという名のオブジェクトとして作業環境内に格納しておこう。
dfは2773行4列のデータであり、2021年9月に実施された世論調査の一部である。IDは回答者の識別番号、Femaleは女性ダミー変数（つまり、女性なら1、それ以外なら0）、Ageは回答者の年齢、Temp_Ishinは日本維新の会に対する感情温度であり、0なら嫌い、100なら好きを意味する。
母平均の推定
まず、母集団（日本人の有権者全体）における維新感情温度の平均値（=母平均）を推定してみよう。母平均の一致推定量&不偏推定量は標本平均（\(\bar{x}\)）である。したがって、Temp_Ishin変数の平均値を計算するだけで、日本人有権者における維新感情温度の平均値が推定できる。ちなみにTemp_Ishinには欠損値が含まれているため、平均値を計算する際、na.rm = TRUEを忘れないこと。また、今後の分析のためにTemp_Ishinの標本分散（不偏分散; \(s^2\)）と有効ケース数（\(n\)）も計算しておこう。計算結果はsum_statという名のオブジェクトとして格納しておく。
sum_stat <- df %>%
summarise(Mean = mean(Temp_Ishin, na.rm = TRUE),
Var = var(Temp_Ishin, na.rm = TRUE),
N = sum(!is.na(Temp_Ishin)))
sum_stat# A tibble: 1 × 3
Mean Var N
<dbl> <dbl> <int>
1 42.7 569. 1998 標本平均は42.733である。これは母集団（日本人の有権者）における維新感情温度の平均値の点推定値が42.733度であることを意味する。しかし、点推定値には常に不確実性が伴う。私たちがもう一回（2021年9月に戻って）同じやり方で2015名を対象とした世論調査を行えば、その時の点推定値（=標本平均）は異なるだろう。この点推定値（=標本平均）の不確実性を示す指標が標準誤差である。標準誤差は\(\sqrt{\frac{s^2}{n}}\)、あるいは\(\frac{s}{\sqrt{n}}\)で計算できる。sum_statに標準誤差をSEという名の列としてN列の前に追加し、sum_statを上書きする。
# A tibble: 1 × 4
Mean Var SE N
<dbl> <dbl> <dbl> <int>
1 42.7 569. 0.534 1998 標準誤差は約0.534であるが、これだけだとどう解釈すれば良いかが分からないかも知れない（ある程度統計学に慣れたら標準誤差だけでも色々な考察ができるようになる）。標準誤差を使って、95%信頼区間を計算してみよう。95%信頼区間を計算するためには「（1）自由度の\(t\)分布において\(T\)統計量が\(t\)以下の値を取る確率が2.5%である\(t\)の値（\(t_{1997, 0.025}\)）」と「（2）自由度の\(t\)分布において\(T\)統計量が\(t\)以上の値を取る確率が2.5%である\(t\)の値（\(t_{1997, 0.975}\)）」を計算する必要がある。使用する確率分布は\(t\)分布であるため、qt()関数を使う。
まず、「（1）自由度の\(t\)分布において\(T\)統計量が\(t\)以下の値を取る確率が2.5%である\(t\)の値」から計算してみよう。qt()の第一引数として0.025を指定する。これは2.5%を意味する。続いて、df引数に\(t\)分布の自由度を指定する。自由度はTemp_Ishinの有効ケース数（\(n\)）から1を引いた値だ。
得られた値は約-1.961である。これは自由度の\(t\)分布において\(T\)統計量が-1.961以下の値を取る確率が2.5%であることを意味する。
つづいて、「（2）自由度の\(t\)分布において\(T\)統計量が\(t\)以上の値を取る確率が2.5%である\(t\)の値」から計算してみよう。これは逆にいうと「\(T\)統計量が\(t\)以上の値を取る確率が97.5%である\(t\)の値」とも解釈できる。したがって、qt()の第一引数として0.975を指定する。
得られた値は約1.961である。これは自由度の\(t\)分布において\(T\)統計量が1.961以上の値を取る確率が2.5%であることを意味する。
95%信頼区間の下限は\(\bar{x} + t_{1997, 0.025} \cdot \mbox{SE}\)、上限は\(\bar{x} + t_{1997, 0.975} \cdot \mbox{SE}\)で計算できる（詳細はスライド参照）。それぞれCI_lower、CI_upperという名の列としてN列前に追加し、sum_statを上書きする。
sum_stat <- sum_stat %>%
mutate(CI_lower = Mean + qt(0.025, df = N - 1) * SE,
CI_upper = Mean + qt(0.975, df = N - 1) * SE,
.before = N)
sum_stat# A tibble: 1 × 6
Mean Var SE CI_lower CI_upper N
<dbl> <dbl> <dbl> <dbl> <dbl> <int>
1 42.7 569. 0.534 41.7 43.8 1998小数点を含めた詳しい値を確認してみよう。
[1] 42.73273[1] 0.5338369[1] 41.6858[1] 43.77967以上の計算により、標本平均は42.733、95%信頼区間は[41.686, 43.780]であることが分かる。
以上の手順はパイプ演算子を複数つなげることで、dfから始まる一つのコードとしてまとめることもできる。複数のパイプ演算子をつなげていくためには、自分の頭の中で処理過程が想像できるようになる必要がある。学習段階では以上のようにステップごとにオブジェクトとして格納し、確認しながら分析を進めることが良いだろう。
# オブジェクトとして格納せず、結果だけを出力
df %>%
summarise(Mean = mean(Temp_Ishin, na.rm = TRUE),
Var = var(Temp_Ishin, na.rm = TRUE),
N = sum(!is.na(Temp_Ishin))) %>%
mutate(SE = sqrt(Var / N),
t = Mean / SE,
CI_lower = Mean + qt(0.025, df = N - 1) * SE,
CI_upper = Mean + qt(0.975, df = N - 1) * SE,
.before = N)# A tibble: 1 × 7
Mean Var SE t CI_lower CI_upper N
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
1 42.7 569. 0.534 80.0 41.7 43.8 1998t.test()の使い方
母平均の推定について学習する場合、これまでのように一つ一つ計算しながら結果を確認する方法が望ましいが、実際にはt.test()という便利な関数を使う。指定する引数は2つであり、第一引数は母平均を推定するベクトル名、第二引数はconf.levelであり、算出する信頼区間を指定する。95%信頼区間の場合、0.95で良い。ちなみにconf.levelの既定値は0.95であるため、95%信頼区間を計算する場合は省略可能である。
One Sample t-test
data: df$Temp_Ishin
t = 80.048, df = 1997, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
41.68580 43.77967
sample estimates:
mean of x
42.73273 手計算（？）と同じ結果が得られることが分かる。