20. 回帰分析がよくわかるシミュレーション
Masahiko Asano
2022-09-10
- このセクションで使っている
R packages
library(tidyverse)
library(stargazer)
library(DT)1. モンテ・カルロシミュレーション
- 回帰分析のしくみ（特に信頼区間の意味）を理解するために、簡単なモンテ・
カルロシミュレーションを行う
- シミュレーションでは、自分で母数を設定し、データを人工的に作り出す
- 特徴を知りたい分析手法（ここでは
OLS）を人工的に作り出したデータに当てはめ、母数を期待通りに推測できるかどうか確認する
- ここでは、単回帰を例にシミュレーションを行う このシミュレーションを行う主な目的は以下の 3 つ
-
最小二乗法が想定するデータ生成過程
(
data generating process) を理解する -
最小二乗推定量（
OLSE: Ordinary Least Squares Estimator）の基本的な性質を理解する -
信頼区間 (
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 の一様分布から長さ
nのxベクトルを無作為抽出
x <- runif(n, 0, 10) yの値を指定
y <- rnorm(n, alpha1 + alpha2 * x, sd = sigma)- 設定に従って正規分布から抽出した
xとyをデータフレームに入れる
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.7597515tail(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 による計量政治学』第７章を参照して下さい）
sd(sim.1$a2)[1] 0.2077991結果 シミュレーションの結果、標準誤差 (se) は推定値の標準偏差 (sd)に（ほぼ）一致する
1.3 信頼区間の理解
- 傾きの係数の推定値 \(α_2\) の
95パーセント信頼区間を例として考える
a2.lowerとa2.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.lowerとa2.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.688604450番目と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 >= alpha2TRUEとなっているものが \(α_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%信頼区間に含まれる結果に違いがあるか確かめなさい
- 総選挙のマップ作成は遠藤勇哉氏（東北大学大学院情報科学研究科博士後期課程）の助言に基づいている
- 宋財泫 (Jaehyun Song)・矢内勇生 (Yuki Yanai)「私たちのR: ベストプラクティスの探究」
- 土井翔平（北海道大学公共政策大学院）「Rで計量政治学入門」
- 矢内勇生（高知工科大学）授業一覧
- 浅野正彦, 矢内勇生.『Rによる計量政治学』オーム社、2018年
- 浅野正彦, 中村公亮.『初めてのRStudio』オーム社、2018年
- Winston Chang, R Graphics Cookbook, O’Reilly Media, 2012.
- Kieran Healy, DATA VISUALIZATION, Princeton, 2019
- Kosuke Imai, Quantitative Social Science: An Introduction, Princeton University Press, 2017