浅野・矢内『Rによる計量政治学』(2018,オーム社)

Chapter 07 統計的推計

浅野・矢内(2018)に基づいて７章の統計的推計をまとめていきます。コマンドは、オーム社のHPからDLしたスクリプトを元にしています。また、他に、石村貞夫『統計解析のはなし』（東京図書、1989年）などを参考にしています。
##7.1　母集団と標本 ###イントロダクション 母集団の数値、母数(parameter)を推定することを目的とする。 母集団から標本抽出し得られたデータを基に、計算して得られるのが統計量(statistic)である。さらにその中でも母数の推定のために利用されるものが推定量(estimator)である。

推定量に望まれる性質

不偏性(unbiasedness)：推定量を平均すると、求めたい母数に一致するという性質。このような推定量は不偏推定量(unbiased estimator)と呼ばれる。
➞標本から得られる標本平均という統計量を使って母平均という母母数を推定する。
一致性(consistency)：サンプルサイズを大きくするほど、求めたい母数から離れた値をとる確率が小さくなるという性質。このような推定量は一致推定量(consistent estimator)と呼ばれる。

7.２　標本分布

まずは今回使う tidyverse パッケージを読み込む

library("tidyverse")
## ─ Attaching packages ───────── tidyverse 1.2.1 ─
## ✔ ggplot2 3.1.0     ✔ purrr   0.3.0
## ✔ tibble  1.4.2     ✔ dplyr   0.7.8
## ✔ tidyr   0.8.2     ✔ stringr 1.3.1
## ✔ readr   1.3.1     ✔ forcats 0.3.0
## Warning: package 'purrr' was built under R version 3.5.2
## ─ Conflicts ────────── tidyverse_conflicts() ─
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()

統計的推計の実験的確認

統計的推計とは、母集団から一部の標本を抜き出し、この標本を分析することで元の母集団について知識を得ようとすることである。
ここでは、その試行として、仮の母集団データを生成し、標本調査を行う。
まず、生成する母集団データは、10万人の身長データである。この母集団の平均（母平均μ）は170cm、母集団の標準偏差（母標準偏差σ）は6とする。μ=170,σ=6

母集団の作成

set.seed(2018 - 11 - 20)  #乱数を準備
population <- rnorm(100000, mean = 170, sd = 6) #rnorm : Normal Disribution(正規分布)を作成するコマンド
mean(population) #平均値
## [1] 170.0057
sd(population)　#標準偏差(Standard Deviation)
## [1] 6.000263

母集団のヒストグラム

次に、母集団の分布をヒストグラムで図示する。(p.125, 図3)

hist_pop <- data_frame( pop = population) %>% 
  ggplot(aes(x = pop, y = ..density..)) + #確率密度(density)オプションを指定。
    geom_histogram(binwidth = 1, color = "black") +
    xlim(140, 200) + ylim(0, 0.6) +
    labs(x = "height (cm)", y = "density")
print(hist_pop) #図示

母集団から標本抽出

上で定義した母集団から、サイズ10の標本を500個抽出する。結果を保存するために、予め10行x500列の空の行列を用意しておく。 新たに使用するのは matrix, for, sample

size10 <- matrix(data = NA, nrow = 10, ncol = 500) #matrixコマンドで、10行×500列の行列を、size10という名前で作成。
for (i in 1:500) {
  size10[, i] <- sample(population, size = 10, replace = FALSE)
# for 構文で同内容の標本抽出を500回繰り返す
#sampleコマンドで、populationというデータフレームの中から、サイズ10のサンプルを、非復元抽出で標本抽出し、size10に格納する。
}

500の標本の平均値

それぞれの標本での平均値を計算し、データフレームの1変数として保存する。
colMeans でsize10の平均値を計算し、means_s10 (10のサンプルサイズの平均値)という名前で、データ保存作成用のdata_frameコマンドで格納。そしてそれをdfというオブジェクト名にする。一応、データを確認。

df <- data_frame(means_s10 = colMeans(size10))
head(df, n = 10) # df 入力で全表示でも良いが、500回の標本抽出の平均値500個になるので長い。
## # A tibble: 10 x 1
##    means_s10
##        <dbl>
##  1      172.
##  2      167.
##  3      168.
##  4      171.
##  5      170.
##  6      172.
##  7      169.
##  8      168.
##  9      172.
## 10      172.

抽出した標本平均のヒストグラム

標本平均の標本分布をヒストグラムにする(p.126 図7.4)

hist_size10 <- ggplot(df, aes(x = means_s10, y = ..density.. )) +
  geom_histogram(binwidth = 1, color = "black") +
  xlim(140, 200) + ylim(0, 0.6) +
  labs(x = expression(paste(" height ", bar(x), " (cm)")),
       y = "density")
print(hist_size10)
## Warning: Removed 2 rows containing missing values (geom_bar).

サンプルサイズ10のデータの確認

500の標本平均内の標準誤差（バラツキ）を計算する

sd(df$means_s10)
## [1] 1.912701

標準偏差: sd(df$means_s10) だとわかる。 ヒストグラムからも標準偏差からも、標本平均が収束していることが伺える。

サイズ90の標本抽出の場合

さらに、サンプルサイズを増やして考えてみる。
先程と同様に、母集団からサンプルサイズ90の標本を500個抽出する。
90行x500列の空の行列を用意し、サンプルサイズ90で標本抽出を500回繰り返す。

size90 <- matrix(NA, nrow = 90, ncol = 500)
for (i in 1:500) { # 標本抽出を500回繰り返す
  size90[, i] <- sample(population, size = 90, replace = FALSE)
}

それぞれの標本での平均値を計算し、データフレームの1変数として保存する

df <- mutate(df, means_s90 = colMeans(size90))
#size90 という名前の標本の平均値を出し、dfの中に、means_s90という名前で変数として新たに作成する。

サンプルサイズ90のヒストグラム

同様に、標本平均の標本分布をヒストグラムにする(p.128,図7.5)

hist_size90 <- ggplot(df, aes(x = means_s90, y = ..density.. )) +
  geom_histogram(binwidth = 1, color = "black") +
  xlim(140, 200) + ylim(0, 0.6) +
  labs(x = expression(paste("height ", bar(x), " (cm)")),
       y = "density")
print(hist_size90)
## Warning: Removed 2 rows containing missing values (geom_bar).

標本サイズ90のデータの確認

標本サイズ90の場合の標本平均のばらつき（標準誤差）を計算する

sd(df$means_s90)
## [1] 0.627559

最後に、確認として、標本サイズ10の場合と90の場合の標本のばらつきの比を計算する

sd(df$means_s10) / sd(df$means_s90)
## [1] 3.047842

7.3 母平均の推定と信頼区間

推定の精度を評価するために、推定の不確実性を明らかにする必要がある。
そこで、信頼区間(confidence interval: CI)を利用する。
信頼区間とは、「母平均μは〇〇以上☓☓以下である」というように、幅をもたせる推定であり、これを利用した推定は区間推定と呼ばれる。

標準誤差

標本統計量の不確実性は標本分布のばらつきとして現れる。 前提として、標本平均の標準誤差の計算式を利用する。サンプルサイズ n と母標準偏差 σ が判明すれば、標本平均の標準誤差がわかる。
しかし、母標準偏差 σ は、未知である。
そこで、まずはこの σ を推定することにする。
そのため、標本平均の標準誤差（standard error:SE）を利用する。そして、これを利用して標本平均を変形させる。そして、 \[ T=\frac{\bar{x}-μ}{SE} \] を得る。このTは、t分布に従っていることが知られれている。さらに、式から標本平均が母平均に一致するときにT=0となる。

t分布

ここで、t分布をRで表示させてみよう。ｔ分布は自由度と呼ばれる数値によって分布の形が変化する。

df_t <- data_frame(x = seq(-3, 3, length.out = 100)) %>% 
  mutate(`1` = dt(x, df = 1),
         `5` = dt(x, df = 5),
         `99` = dt(x, df = 99)) %>% 
  gather(key = "df", value = "density", `1`:`99`)
t_dist <- ggplot(df_t, aes(x = x, y = density,
                           color = df, linetype = df)) +
  geom_line() +
  labs(x = "", y = "density") +
  scale_x_continuous(breaks = -3:3) +
  scale_color_discrete(name = "degree of freedom") + 
  scale_linetype_discrete(name = "degree of freedom") +
  guides(shape = guide_legend(reverse = TRUE)) 
print(t_dist)

t分布を利用した推定

標本平均の推定では、自由度がn-1のt分布を利用する。
それでは、自由度99のt分布で、データの95%が収まる範囲の上限値を求める

qt(p = 0.025, df = 99, lower.tail = FALSE)
## [1] 1.984217
# または、
qt(p = 0.975, df = 99, lower.tail = TRUE)
## [1] 1.984217
## 同様に、下限値を求める
qt(p = 0.025, df = 99, lower.tail = TRUE)
## [1] -1.984217

ここから、データの95％が収まる範囲が、-.198~1.98であることがわかる。

n=100の身長データで試行

上で定義した母集団（身長の母平均が約170）からサンプルサイズ100の標本を1000個抽出し、それぞれの標本について信頼区間を求める。

n <- 100
res <- matrix(NA, nrow = n, ncol = 1000)
for (i in 1:1000) {
  res[, i] <- sample(population, size = n, replace = FALSE)
}

それぞれの標本の平均と標準偏差を求めてデータフレームの変数にする。
その後、標準誤差 se と95%信頼区間の下限 lb、上限ub を計算する。

df_ci <- data_frame(id = 1:1000,
                    means = colMeans(res),
                    sds = apply(res, MARGIN = 2, FUN = sd)) %>% 
  mutate(se = sds / sqrt(n),
         lb = means + qt(p = 0.025, df = 99, lower.tail = TRUE) * se,
         ub = means + qt(p = 0.025, df = 99, lower.tail = FALSE) * se)
head(df_ci) # 結果の一部を見てみる
## # A tibble: 6 x 6
##      id means   sds    se    lb    ub
##   <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1     1  170.  5.40 0.540  169.  171.
## 2     2  170.  5.59 0.559  168.  171.
## 3     3  169.  6.64 0.664  168.  171.
## 4     4  170.  5.69 0.569  169.  171.
## 5     5  169.  6.03 0.603  168.  170.
## 6     6  171.  6.42 0.642  169.  172.

95%信頼区間を求める

信頼区間の中に真の母平均が含まれているかどうか判定する変数を作る

(pop_mean <- mean(population)) # 真の母平均をpop_meanで格納
## [1] 170.0057
df_ci <- mutate(df_ci,
                success = (lb <= pop_mean & ub >= pop_mean))
## 母平均を区間内に捉えた95%信頼区間の割合を計算する
mean(df_ci$success)
## [1] 0.949

オマケ

1000個の標本から50個をランダムに選び、95%信頼区間を図示する。

set.seed(2018-11-22)
ci_plt <- sample_n(df_ci, size = 50) %>% 
  ggplot(aes(x = reorder(id, means), y = means, 
             ymin = lb, ymax = ub, color = success)) +
  geom_hline(yintercept = pop_mean, color = "dodgerblue") +
  geom_linerange() +
  geom_point() +
  scale_color_discrete(name = "does this include the population mean??") +
  labs(x = "", y = "95% confidence interval") +
  coord_flip()
print(ci_plt)

衆院選データでもやってみよう

例7.4：衆院データを使って信頼区間を求める

HR <- read_rds("data/hr-data.Rds")
## Rdsファイルの読み込みがうまくいかない場合は以下を実行
#download.file(url = "https://git.io/fxhQU",
#              destfile = "data/hr-data.csv")
#HR <- read_csv("data/hr-data.csv")
## 候補者の年齢ageの50%信頼区間を求める
age_ci50 <- t.test(HR$age, conf.level = 0.5)
age_ci50$conf.int
## [1] 50.81778 50.97717
## attr(,"conf.level")
## [1] 0.5
## 候補者の年齢ageの95%信頼区間を求める
age_ci95 <- t.test(HR$age)
age_ci95$conf.int
## [1] 50.66587 51.12908
## attr(,"conf.level")
## [1] 0.95

50％信頼区間よりも、95%信頼区間のほうが長い範囲であることがわかる。