第5回講義資料

統計的推定

スライド

セットアップ

　今回使用するパッケージは{tidyverse}のみだ。データはinference_data.csvであり、LMSから入手できる。読み込んだデータはdfという名のオブジェクトとして作業環境内に格納しておこう。

Code 01

library(tidyverse) # or pacman::p_load(tidyverse)
df <- read_csv("Data/inference_data.csv")

df

# A tibble: 2,746 × 4
      ID Female   Age Temp_Ishin
   <dbl>  <dbl> <dbl>      <dbl>
 1     1      0    44         NA
 2     2      1    32         50
 3     3      1    53         20
 4     4      1    22          0
 5     5      1    27         NA
 6     6      1    28         NA
 7     7      0    28         30
 8     8      1    45         NA
 9     9      1    33         50
10    10      0    30         NA
# ℹ 2,736 more rows

　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という名のオブジェクトとして格納しておく。

Code 02

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を上書きする。

Code 03

sum_stat <- sum_stat |>
  mutate(SE      = sqrt(Var / N),
         .before = 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を引いた値だ。

Code 04

qt(0.025, df = sum_stat$N - 1)

[1] -1.961153

　得られた値は約-1.961である。これは自由度の\(t\)分布において\(T\)統計量が-1.961以下の値を取る確率が2.5%であることを意味する。

　つづいて、「（2）自由度の\(t\)分布において\(T\)統計量が\(t\)以上の値を取る確率が2.5%である\(t\)の値」から計算してみよう。これは逆にいうと「\(T\)統計量が\(t\)以上の値を取る確率が97.5%である\(t\)の値」とも解釈できる。したがって、qt()の第一引数として0.975を指定する。

Code 05

qt(0.975, df = sum_stat$N - 1)

[1] 1.961153

　得られた値は約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を上書きする。

Code 06

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

　小数点を含めた詳しい値を確認してみよう。

Code 07

sum_stat$Mean # 母平均の一致＆不偏推定量（標本平均）

[1] 42.73273

Code 08

sum_stat$SE # 標準誤差（標本平均の不確実性）

[1] 0.5338369

Code 09

sum_stat$CI_lower # 95%信頼区間の下限

[1] 41.6858

Code 10

sum_stat$CI_upper # 95%信頼区間の上限

[1] 43.77967

　以上の計算により、標本平均は42.733、95%信頼区間は[41.686, 43.780]であることが分かる。

　以上の手順はパイプ演算子を複数つなげることで、dfから始まる一つのコードとしてまとめることもできる。複数のパイプ演算子をつなげていくためには、自分の頭の中で処理過程が想像できるようになる必要がある。学習段階では以上のようにステップごとにオブジェクトとして格納し、確認しながら分析を進めることが良いだろう。

Code 11

# オブジェクトとして格納せず、結果だけを出力
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  1998

`t.test()`の使い方

　母平均の推定について学習する場合、これまでのように一つ一つ計算しながら結果を確認する方法が望ましいが、実際にはt.test()という便利な関数を使う。指定する引数は2つであり、第一引数は母平均を推定するベクトル名、第二引数はconf.levelであり、算出する信頼区間を指定する。95%信頼区間の場合、0.95で良い。ちなみにconf.levelの既定値は0.95であるため、95%信頼区間を計算する場合は省略可能である。

Code 12

# 95%信頼区間を出すのであれば、 t.test(df$Temp_Ishin) だけでもOK
t.test(df$Temp_Ishin, conf.level = 0.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

　手計算（？）と同じ結果が得られることが分かる。

スライド

セットアップ

母平均の推定

t.test()の使い方

`t.test()`の使い方