• このセクションで使っている R packages
library(tidyverse)
library(stargazer)
library(DT)

1. モンテ・カルロシミュレーション

  • 回帰分析のしくみ(特に信頼区間の意味)を理解するために、簡単なモンテ・ カルロシミュレーションを行う 
  • シミュレーションでは、自分で母数を設定し、データを人工的に作り出す 
  • 特徴を知りたい分析手法(ここでは OLS )を人工的に作り出したデータに当てはめ、母数を期待通りに推測できるかどうか確認する 
  • ここでは、単回帰を例にシミュレーションを行う  このシミュレーションを行う主な目的は以下の 3 つ 
  1. 最小二乗法が想定するデータ生成過程 (data generating process) を理解する
  2. 最小二乗推定量(OLSE: Ordinary Least Squares Estimator)の基本的な性質を理解する
  3. 信頼区間 (Confidence Interval) の意味を理解する
  • 単回帰モデルは次の式で表せる  

\[y_i \sim \mathrm{N}(\alpha_1 + \alpha_2 x_i, \sigma^2)- - - i = 1, 2, \dots, n\]

  • 設定する母数は次の 3 つである 
alpha1 <- 10 # 切片の値は 10  
alpha2 <- 3  # 傾き(独立変数の係数)は 3  
sigma  <- 6  # 標準偏差は 6  
  • 単回帰モデルが想定するデータ生成過程*に従って、データを生成する 

  • サンプルサイズ = 100 と指定

n <- 100   
  • 同一のサンプル結果を得られるよう特定のノイズを設定
set.seed(20210328)
  • 最大値 10 の一様分布から長さ nx ベクトルを無作為抽出
x <- runif(n, 0, 10) 
  • y の値を指定
y <- rnorm(n, alpha1 + alpha2 * x, sd = sigma)
  • 設定に従って正規分布から抽出した xy をデータフレームに入れる
df <- data.frame(y = y, x = x)
head(df)   # データの文頭を表示
          y         x
1 31.797051 8.1281966
2  8.285991 1.7148744
3 20.963492 1.9013890
4 29.246675 6.4274074
5 33.201243 9.5686632
6 12.079381 0.7597515
tail(df)  # データの終わりを表示
           y        x
95  23.88150 8.349453
96  25.52376 6.728700
97  42.03960 9.139047
98  34.03598 7.244872
99  25.14174 5.804608
100 16.32483 3.505625
  • 散布図にして直線を当てはめる
library("ggplot2")
scat <- ggplot(df, aes(x, y)) + 
  geom_point() + 
  geom_smooth(method = "lm")
  print(scat)

  • 最初に母数として次の 3 つの値を設定した: 

\[切片:α_1=10\] \[傾き:α_2=3\] \[標準偏差:sigma=6\]

  • このサンプルから推定値を計算する:
reg1 <- lm(y ~ x, data = df)
summary(reg1)

Call:
lm(formula = y ~ x, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-19.5307  -3.4985   0.3789   4.2163  11.0789 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  11.5055     1.0480   10.98   <2e-16 ***
x             2.7065     0.1852   14.62   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.737 on 98 degrees of freedom
Multiple R-squared:  0.6856,    Adjusted R-squared:  0.6823 
F-statistic: 213.7 on 1 and 98 DF,  p-value: < 2.2e-16

\[切片: Intercept = 11.5055 \\ 傾き: Coefficients = 2.7065 \\ 標準誤差: Standard.error = 5.737\]

  • それぞれの母数 (10, 3, 6) に近い値が得られた 

  • confint(reg1) とコマンドを入力すると係数の95パーセント信頼区間が得られる 

confint(reg1)
               2.5 %   97.5 %
(Intercept) 9.425757 13.58534
x           2.339090  3.07399
  • この 95 パーセント信頼区間は最初に設定した母数 \(α_1=10, α_2=3\) を含んでいる 

ここで得られた信頼区間は 100% の確率で母数を含む

  • 以上の過程を、母数とサンプルサイズは変えずに何度も繰り返す 
  • 同じコードを何度も繰り返し使うのは面倒なので、関数にする 
sim_ols1 <- function(alpha, sigma,   # α は傾き、σ は標準偏差
                     n = 100,     # サンプルサイズ
                     trials = 10000, # 繰り返し回数を指定 
                     x.rng = c(0, 10)) { # 説明変数 x の範囲指定
    

    empty <- rep(NA, trials)
    df <- data.frame(   # df に結果を保存
        a1 = empty,   # データフレームには次の 9 つの変数を含める
        a1.se = empty,
        a1.lower = empty,
        a1.upper = empty,
        a2 = empty,
        a2.se = empty,
        a2.lower = empty,
        a2.upper = empty,
        sigma.hat = empty
        )
    
    for (i in 1:trials) { ## loop for trials
        ## generate data
        x <- runif(n, x.rng[1], x.rng[2])
        y <- rnorm(n, alpha[1] + alpha[2] * x, sigma)
        
        ## fit OLS
        fit <- lm(y ~ x)
        
        ## sum of squared residuals
        sigma.hat <- summary(fit)$sigma
        
        ## coefficient estimates
        a1 <- coef(fit)[1]
        a2 <- coef(fit)[2]
        
        ## standard errors of coefficient estimates
        a1.se <- sqrt(summary(fit)$cov.unscaled[1,1]) * sigma.hat
        a2.se <- sqrt(summary(fit)$cov.unscaled[2,2]) * sigma.hat
        
        ## 95% CI
        a1.ci95 <- confint(fit)[1,]
        a2.ci95 <- confint(fit)[2,]
        
        ## save the results
        df[i,] <- c(a1, a1.se, a1.ci95,
                    a2, a2.se, a2.ci95,
                    sigma.hat)      
    }
    ## return the data frame
    return(df)
}
  • この関数を利用して、実際にシミュレーションを行う 
  • 母数とサンプルサイズ(ここでは規定値の n = 100 を利用)を指定し、データの生成と回帰式の推定を1,000回繰り返してみる 
alpha1 <- 10; 
alpha2 <- 3; 
sigma = 6
sim.1 <- sim_ols1(alpha = c(alpha1, alpha2), 
                  sigma = sigma, 
                  trials = 1000)
  • このシミュレーションで得られた結果は次のとおり
summary(sim.1)
       a1             a1.se           a1.lower         a1.upper     
 Min.   : 5.640   Min.   :0.8696   Min.   : 3.282   Min.   : 7.613  
 1st Qu.: 9.152   1st Qu.:1.1216   1st Qu.: 6.738   1st Qu.:11.527  
 Median : 9.890   Median :1.1940   Median : 7.532   Median :12.314  
 Mean   : 9.946   Mean   :1.1996   Mean   : 7.566   Mean   :12.327  
 3rd Qu.:10.772   3rd Qu.:1.2726   3rd Qu.: 8.398   3rd Qu.:13.177  
 Max.   :15.434   Max.   :1.5432   Max.   :12.900   Max.   :17.968  
       a2            a2.se           a2.lower        a2.upper    
 Min.   :2.255   Min.   :0.1579   Min.   :1.814   Min.   :2.696  
 1st Qu.:2.863   1st Qu.:0.1949   1st Qu.:2.448   1st Qu.:3.268  
 Median :3.006   Median :0.2081   Median :2.590   Median :3.417  
 Mean   :3.008   Mean   :0.2080   Mean   :2.595   Mean   :3.421  
 3rd Qu.:3.146   3rd Qu.:0.2201   3rd Qu.:2.737   3rd Qu.:3.565  
 Max.   :3.725   Max.   :0.2739   Max.   :3.405   Max.   :4.082  
   sigma.hat    
 Min.   :4.619  
 1st Qu.:5.645  
 Median :5.953  
 Mean   :5.962  
 3rd Qu.:6.257  
 Max.   :7.190  
変数名 詳細
a1 切片 (Intercept) の推定値
a1.se 切片 (Intercept) の標準誤差
a1.lower 切片の95パーセント信頼区間(下側)
a1.upper 切片の95パーセント信頼区間(上側)
a2 傾き (Coefficients) の推定値
a2.se 傾き (Coefficients)の標準誤差
a2.lower 傾きの95パーセント信頼区間(下側)
a2.upper 傾きの95パーセント信頼区間(上側)
sigma.hat 標準偏差 (standard error) の推定値
DT::datatable(sim.1)

1.1 偏回帰係数(傾き)の理解

  • 上で得れたシミュレーションの結果 sim.1 を使って、説明変数 x の係数の推定値 \(α_2\) (傾き)の分布を確認してみよう 
  • 私たちは母数の傾き \(α_2 = 3\)であることを知っているが、OLS はこの数値をうまく推定してくれただろうか?
hist.a2 <- ggplot(sim.1, 
                  aes(x = a2,
                      y =..density..)) + 
  geom_histogram(binwidth = .1, color = "white")

hist.a2 <- hist.a2 + labs(x = "傾き",
                          y = "確率密度", 
           title = "シミュレーション結果:傾き(α2)の分布") +  
           geom_vline(xintercept = 3, color = "green") + # 傾き= 3 に縦線を引く
  theme_bw(base_family = "HiraKakuProN-W3") 

print(hist.a2)

  • 上記のヒストグラムが示すように、OLS は平均すると、母数として設定した \(α_2 = 3\) をうまく推定している 
  • これを不偏性 (unbiasednessという 
  • しかし、常に正しい推定をするわけではない 
  • 係数の値を過小推定することもあれば、過大推定することもある 
  • 実際の分析では、データセットが1つしかないのが普通 
     →自分のデータ分析が係数を「正しく」推定しているとは限らならい 
  • そのために、推定の不確実性を明示することが求められる 

1.2 標準誤差の理解

  • \(α_2\) の標準誤差をヒストグラムで表示する 
hist.se <- ggplot(sim.1, 
                  aes(x = a2.se,      # x 軸に表示するのは α2 の標準誤差
                      y =..density..)) + 
    geom_histogram(binwidth = .01, 
                   color = "white") +
    labs(x = "alpha_2(傾き)の標準誤差", 
         y = "確率密度", 
         title = "シミュレーション結果:alpha_2(傾き)の標準誤差") +
  theme_bw(base_family = "HiraKakuProN-W3") 
print(hist.se)

  • 標準誤差 (se) 自体が推定量なので、その値はサンプルごとに異なる(つまり、分布する)   

  • 標準誤差 (se) は次の式で求められる 

\[\frac{\hat{\alpha} - \alpha}{\mathrm{se}}\]

  • これは自由度 \((n-k)\)(ここでは、100 - 2 = 98)の \(t\) 分布に従うはず    
  • ここで求めた標準誤差 (se) をヒストグラムにして、自由度98の \(t\) 分布の確率密度曲線を重ね書きするしてみる 
library("dplyr")
sim.1 <- mutate(sim.1, a2.t = (a2 - alpha2) / a2.se)
    ## t分布の確率密度を求め、true.t という名前のデータフレームに保存する
true.t <- data.frame(x = seq(-4, 4, length = 100))
true.t <- mutate(true.t, density = dt(x, df = 98))
hist.t <- ggplot(sim.1, aes(x = a2.t, y =..density..)) +
    geom_histogram(binwidth = 0.5, color = "white") + 
    geom_line(data = true.t, aes(x = x, y = density), color = "violetred") +
  theme_bw(base_family = "HiraKakuProN-W3") 

hist.t <- hist.t +
    labs(x = expression(frac(hat(alpha) - alpha, "se")),
         y="確率密度", title="標準化した α2 と t 分布 (df = 98)")
print(hist.t)

  • 上記の図から、標準誤差 (se) の分布は \(t\) 分布に近いことがわかる 
  • Q-Q プロットでも確かめてみよう 
qqplot.t <- ggplot(sim.1, aes(sample = a2.t)) + 
    stat_qq(distribution = qt, dparams = list(df = 98)) +
    geom_abline(intercept = 0, slope = 1, color = "violetred") +
    labs(title = "Q-Q Plot", 
         x = "t 分布 (df = 98)",
         y = expression(paste("シミュレーション結果 ", frac(hat(alpha) - alpha, "se"), ))) +
  theme_bw(base_family = "HiraKakuProN-W3") 
print(qqplot.t)

  • Q-Q プロットの点がほぼ一直線に並んでいる 

  • \((\hat{\alpha}-\alpha)/\mathrm{se})\)\(t\) 分布に従っている 

  • ただし、分布の裾では、理論値との乖離が大きいことに注意 

  • 今回のシミュレーションで得られた標準誤差の平均値をもとめてみる   

mean(sim.1$a2.se)
[1] 0.2080031
  • 標準誤差 (se) は推定値 \(α_1\) の標準偏差のはずなので (標準偏差と標準誤差の関係については『R による計量政治学』第7章を参照して下さい)
sd(sim.1$a2)
[1] 0.2077991

結果 シミュレーションの結果、標準誤差 (se) は推定値の標準偏差 (sd)に(ほぼ)一致する

1.3 信頼区間の理解

  • 傾きの係数の推定値 \(α_2\) の 95パーセント信頼区間を例として考える 
  • a2.lowera2.upper の間が \(α_2\) の95パーセント信頼区間 
  • 上記のシミュレーションから傾きの係数の推定値 \(α_2\) の「95パーセント信頼区間」に関しては次の結果が得られた:
  • 結果の最初の 3 行だけ表示してみる
sim.1[1:3, c("a2.lower", "a2.upper")]
  a2.lower a2.upper
1 2.651229 3.528330
2 2.415444 3.117649
3 2.482895 3.305901
  • 結果の最後の 3 行だけ表示してみる
sim.1[998:1000, c("a2.lower", "a2.upper")]
     a2.lower a2.upper
998  2.659470 3.444083
999  2.553984 3.287236
1000 2.746300 3.519873
  • 母数である \(α_2 = 3\)a2.lowera2.upper の区間に含まれる 
  • 例えば、450番目 から452番目までの信頼区間を表示する 
sim.1[450:452, c("a2.lower", "a2.upper")]
    a2.lower a2.upper
450 2.367338 3.257226
451 2.016128 2.887802
452 2.914034 3.688604
  • 450番目と452番目は母数 \(α_2 = 3\) を区間内に含んでいる
    →この信頼区間が母数 \(α_2 = 3\) を含む確率は 100% 

  • しかし451番目は母数 \(α_2 = 3\) を区間内に含んでいない
    →この信頼区間が母数 \(α_2 = 3\) を含む確率は 0%

  • シミュレーションで得た 1,000 個の信頼区間のうち、どれだけの信頼区間が母数母数 \(α_2 = 3\) を含んでいるか調べてみる 

  • 母数 \(α_2 = 3\) を信頼区間内に含むのは、信頼区間の下限値が母数以下かつ上限値が母数以上のもの 

check.ci <- sim.1$a2.lower <= alpha1 & sim.1$a2.upper >= alpha2
  • TRUE となっているものが \(α_2 = 3\) を区間内に捉えているもの、FALSE がそうでないもの 
  • これを表にすると次のようになる 
table(check.ci)
check.ci
FALSE  TRUE 
   23   977 
  • 1000個の 95 パーセント信頼区間のうち、971個は母数 \(α_2 = 3\) を区間内に捉えたが、29個は捉えていない 
  • 一定のデータ生成過程から得られた異なるデートセットに対し、信頼区間を求める作業を繰り返した時「得られた信頼区間のうち何パーセントが母数を区間内に含むか」というのが「信頼区間の信頼度」
     → 例えば「5% 以内含む」のであれば「5%で統計的に有意」という
  • この結果を図にしてみよう 
  • 1000個は多すぎるので、無作為に 30 個だけ選んで図示する 
sim.1 <- mutate(sim.1, id = 1:n()) # シミュレーションごとに ID 番号を設定
sim.1$check.ci <- check.ci          # check.ci をデータフレームに入れる  
sim.1.sml <- sample_n(sim.1, 30)   # ランダムに 30個の結果を選ぶ

ct_plot <- ggplot(sim.1.sml, 
                aes(x = reorder(id, a2), 
                    y = a2, 
                    ymin = a2.lower,
                    ymax = a2.upper, 
                    color = check.ci)) +
  geom_hline(yintercept = alpha2, 
             linetype = "dashed") + 
  geom_pointrange(size = 0.2) + 
  labs(x = "シミュレーション ID", 
       y = "a2(傾き)",
       title = "95% 信頼区間") +
  scale_color_discrete(name = "mu = 3", 
                       labels = c("95%に含まれない","95%に含まれる")) +
  coord_flip()+
  theme_bw(base_family = "HiraKakuProN-W3") 

print(ct_plot)

  • 30回のうち 2 回 (586番目と451番目) が95%信頼区間を捉えることができず「95%に含まれない」ことがわかる
  • さらに詳しい説明は『R による計量政治学』第7章「統計的推定」を参照して下さい

2. Exercise

  • ここでは同一のサンプル結果を得られるよう特定のノイズを set.seed(20210328)と設定したが、この設定をせず、シミュレーションを実行して可視化し、実行するたびにどの程度母数が95%信頼区間に含まれる結果に違いがあるか確かめなさい
参考文献