計量経済学応用
重回帰分析 (1)
矢内 勇生
2018-04-25 (rev: 2018-04-26)0- 2
準備
まず、必要なパッケージを読み込む。
library("foreign")
library("dplyr")
library("arm")
library("ggplot2")
# Macユーザは次の行も実行する
#theme_set(theme_gray(base_size = 12, base_family = "HiraKakuProN-W3"))前回同様、『Stataによる計量政治学』（浅野正彦, 矢内勇生. 2013）で使用されているデータ（hr96-09.dta） を使う。
## year ku kun party name age status nocand wl rank
## 1 1996 aichi 1 NFP KAWAMURA, TAKASHI 47 incumbent 7 win 1
## 2 1996 aichi 1 LDP IMAEDA, NORIO 72 moto 7 lose 2
## 3 1996 aichi 1 DPJ SATO, TAISUKE 53 incumbent 7 lose 3
## 4 1996 aichi 1 JCP IWANAKA, MIHOKO 43 challenger 7 lose 4
## 5 1996 aichi 1 bunka ITO, MASAKO 51 challenger 7 lose 5
## 6 1996 aichi 1 NP YAMADA, HIROSHIB 51 challenger 7 lose 6
## previous vote voteshare eligible turnout exp
## 1 2 66876 40.0 346774 49.22 9828097
## 2 3 42969 25.7 346774 49.22 9311555
## 3 2 33503 20.1 346774 49.22 9231284
## 4 0 22209 13.3 346774 49.22 2177203
## 5 0 616 0.4 346774 49.22 NA
## 6 0 566 0.3 346774 49.22 NA## year ku kun party name age status nocand
## 5609 2009 yamanashi 2 msz NAGASAKI, KOTARO 41 incumbent 4
## 5610 2009 yamanashi 2 LDP HORIUCHI, MITSUO 79 incumbent 4
## 5611 2009 yamanashi 2 H MIYAMATSU, HIROYUKI 69 challenger 4
## 5612 2009 yamanashi 3 DPJ GOTO, HITOSHI 52 incumbent 3
## 5613 2009 yamanashi 3 LDP ONO, JIRO 56 incumbent 3
## 5614 2009 yamanashi 3 H SAKURADA, DAISUKE 47 challenger 3
## wl rank previous vote voteshare eligible turnout exp
## 5609 lose 2 1 57213 32.1 234746 77.09 7916556
## 5610 lose 3 10 52773 29.6 234746 77.09 11611677
## 5611 lose 4 0 1214 0.7 234746 77.09 1326378
## 5612 win 1 3 112894 62.7 248102 74.70 6795969
## 5613 lose 2 1 63611 35.3 248102 74.70 12876644
## 5614 lose 3 0 3663 2.0 248102 74.70 1953819衆議院議員経験があることを表すダミー変数と選挙費用を100万円単位で測定する変数を作る。
HR <- HR %>%
mutate(experience = as.numeric(status == "incumbent" | status == "moto"),
expm = exp / 10^6)2009年の結果だけ抜き出し、HR09として保存する。
Rで重回帰分析を行う
複数の説明変数を使った回帰（モデル3）
得票率を議員経験と選挙費用で説明するモデルを推定する。 lm() で重回帰を行うときは、説明変数を + でつなぐ。
## lm(formula = voteshare ~ experience + expm, data = HR09)
## coef.est coef.se
## (Intercept) 7.91 0.69
## experience 18.10 1.23
## expm 1.85 0.12
## ---
## n = 1124, k = 3
## residual sd = 14.69, R-Squared = 0.56この結果を図示する。
## 予測値のデータフレームを作る
pred <- with(HR09, expand.grid(
expm = seq(min(expm, na.rm = TRUE), max(expm, na.rm = TRUE), length.out = 100),
experience = c(0,1)
))
## 予測値を計算する
pred <- mutate(pred, voteshare = predict(fit_3, newdata = pred))
## 散布図 (point) は元データ (HR09) で描き、回帰直線は予測値 (pred) で描く
p3 <- ggplot(HR09, aes(x = expm, y = voteshare, color = as.factor(experience))) +
geom_point(size = 1) +
geom_line(data = pred)
p3 <- p3 + labs(x = '選挙費用（100万円）', y = '得票率（%）')
p3 <- p3 +
scale_color_discrete(name = '議員経験', labels = c('なし', 'あり')) +
guides(color = guide_legend(reverse = TRUE))
print(p3 + ggtitle('得票率と選挙費用の関係'))
この図からわかるように、このモデルは2つの直線が平行になる（傾きが同じになる）ように設定されている。
この図に95パーセント信頼区間を加えよう。 信頼区間を求めるために標準誤差を利用する。 また、標準誤差は\(t\) 分布に従うので、qt() で分布の95%が収まる範囲を求める （標準正規分布で近似し、\(\pm 1.96\) を使っても結果に大きな違いはないが、せっかくRを使っているのだから、より正確な数値を利用したほうがよい）。
## 信頼区間の下限と上限を計算する
## まず、予測値と標準誤差を求める
err <- predict(fit_3, newdata = pred, se.fit = TRUE)
## 予測値と標準誤差を使って信頼区間を求める
pred$lower <- err$fit + qt(0.025, df = err$df) * err$se.fit
pred$upper <- err$fit + qt(0.975, df = err$df) * err$se.fit
p3_ci95 <- p3 +
geom_smooth(data = pred, aes(ymin = lower, ymax = upper), stat = 'identity')
print(p3_ci95 + ggtitle('得票率と選挙費用の関係'))
相互作用を考える（モデル4）
モデル3では回帰直線が平行になるようにモデルを制限していたが、この制限を緩めてみよう。 そのために、議員経験 \(\times\) 選挙費用という説明変数を追加する。
Rでは、lm() 内の説明変数のところに、変数A * 変数B と書くと、変数A、変数B、変数A \(\times\) 変数B を説明変数として加えてくれるので、これを利用する。
## lm(formula = voteshare ~ experience * expm, data = HR09)
## coef.est coef.se
## (Intercept) -2.07 0.72
## experience 45.91 1.58
## expm 4.87 0.16
## experience:expm -4.76 0.21
## ---
## n = 1124, k = 4
## residual sd = 12.11, R-Squared = 0.70この結果を図示する。
## 議員経験の有無によって、色分けを行う
p4 <- ggplot(HR09, aes(x = expm, y = voteshare, color = as.factor(experience))) +
geom_point(size = 1) +
geom_smooth(method = 'lm', se = FALSE)
p4 <- p4 + labs(x = '選挙費用（100万円）', y = '得票率（%）')
p4 <- p4 +
scale_color_discrete(name = '議員経験', labels = c('なし', 'あり')) +
guides(color = guide_legend(reverse = TRUE))
print(p4 + ggtitle('得票率と選挙費用の関係'))
95パーセント信頼区間を加えよう。
シミュレーションで回帰分析を理解する
回帰分析のしくみ（特に信頼区間の意味）を理解するために、簡単なモンテカルロシミュレーションを行おう。 シミュレーションでは、自分で母数を設定し、データを生成する。 そのうえで、特徴を知りたい分析手法（今回の場合はOLS）を生成したデータに当てはめ、母数をうまく推測できるかどうか確認する。
今回は、単回帰を例にシミュレーションを行う。 このシミュレーションを行う主な目的は以下の3つである。
- 線形回帰が想定するデータ生成過程 (data generating process) を理解する
- 線形回帰の推定量の基本的な性質を理解する
- 信頼区間の意味を理解する
単回帰モデルは、以下のとおりである。 \[y_i = \beta_1 + \beta_2 x_i + u_i,\] \[u_i \sim \mathrm{N}(0, \sigma^2).\] ここで、\(u_i\) は誤差項である。
あるいは、以下のように表しても同じ意味である。 \[y \sim \mathrm{N}(X\beta, \sigma^2 I)\] または、 \[y_i \sim \mathrm{N}(\beta_1 + \beta_2 x_i, \sigma^2)\] (\(i = 1, 2, \dots, n\))と表すことができる。
したがって、設定する母数は3つ（\(\beta_1\), \(\beta_2\), \(\sigma\)）ある。
次に、単回帰モデルが想定するデータ生成過程に従って、データを生成する。
n <- 100 ## サンプルサイズを決める
## x ~ U(0, 10) とする：他の設定でもかまわない
x <- runif(n, min = 0, max = 10) ## 最小値0, 最大値10の一様分布から長さnのxベクトルを抽出する
y <- rnorm(n, beta1 + beta2 * x, sd = sigma) ## 設定に従って、yを正規分布から抽出するここで、手に入れたデータを散布図にして直線を当てはめてみよう。
df <- data.frame(y, x)
scat <- ggplot(df, aes(x, y)) + geom_point() + geom_smooth(method = 'lm')
print(scat)
回帰式を推定すると、
## lm(formula = y ~ x, data = df)
## coef.est coef.se
## (Intercept) 9.14 1.28
## x 3.10 0.22
## ---
## n = 100, k = 2
## residual sd = 6.29, R-Squared = 0.683つの母数\(\beta_1=10, \beta_2=3, \sigma=6\) に対する推定値は、それぞれ9.14, 3.1, 6.29 であることがわかる。
このとき、係数の95パーセント信頼区間は、
## 2.5 % 97.5 %
## (Intercept) 6.602367 11.686254
## x 2.672670 3.526034で求められる。 95パーセント信頼区間は、切片が[6.60, 11.69]、傾きが[2.67, 3.53] であり、どちらの信頼区間も母数を含んでいる。 つまり、ここで得られた信頼区間が母数を含む確率は1（100%）である！
以上の過程を、母数とサンプルサイズは変えずに何度も繰り返す。 同じコードを何度も繰り返し使うのは面倒なので、関数にしてしまおう。
sim_ols1 <- function(beta, sigma, n = 100, trials = 10000, x_rng = c(0,10)) {
## Function to simulate simple linear regressions
## Arguments:
## beta: vector of coefficient parameters
## sigma: standard deviation of the error
## n: sample size, default n=100
## trials: the number of simulation trials
## default trials=10000
## x_rng: range of the explanatory variable x
## default x.rng=c(0,10)
## Return:
## df: data frame containing
## (1) parameter estimates,
## (2) se of coefficient estimates, and
## (3) 95% CI of coefficient estimates
## for each trial of the simulation.
## create a data frame to save the results
col_names <- c('b1', 'b1_se', 'b1_lower', 'b1_upper',
'b2', 'b2_se', 'b2_lower', 'b2_upper', 'sigma_hat')
df <- as.data.frame(matrix(rep(NA, trials * length(col_names)),
ncol = length(col_names)))
names(df) <- col_names
for (i in 1:trials) { # loop for trials
# generate data
x <- runif(n, x_rng[1], x_rng[2])
y <- rnorm(n, mean = beta[1] + beta[2]*x, sd = sigma)
# regress y on x
fit <- lm(y ~ x)
# sum of squared residuals
sigma_hat <- summary(fit)$sigma
# coefficient estimates
b1 <- coef(fit)[1]
b2 <- coef(fit)[2]
# standard errors of coefficient estimates
b1_se <- sqrt(summary(fit)$cov.unscaled[1,1]) * sigma_hat
b2_se <- sqrt(summary(fit)$cov.unscaled[2,2]) * sigma_hat
# 95% CI
b1_ci95 <- confint(fit)[1,]
b2_ci95 <- confint(fit)[2,]
# save the results
df[i,] <- c(b1, b1_se, b1_ci95,
b2, b2_se, b2_ci95,
sigma_hat)
}
# return the data frame
return(df)
}この関数を利用して、実際にシミュレーションを行ってみよう。 自分で母数とサンプルサイズ（ここでは規定値のn=100を利用）を指定し、データの生成と回帰式の推定を1,000回繰り返すことにする。
beta1 <- 10; beta2 <- 3; sigma <- 6
sim_1 <- sim_ols1(beta = c(beta1, beta2), sigma = sigma, trials = 1000)
head(sim_1)## b1 b1_se b1_lower b1_upper b2 b2_se b2_lower
## 1 10.519381 1.333821 7.872455 13.16631 2.840636 0.2276824 2.388807
## 2 9.195160 1.558476 6.102415 12.28790 3.157845 0.2427141 2.676187
## 3 8.675866 1.365326 5.966421 11.38531 3.169810 0.2311890 2.711023
## 4 9.693470 1.422950 6.869673 12.51727 2.911867 0.2347170 2.446079
## 5 8.809144 1.111461 6.603487 11.01480 3.040280 0.2021615 2.639097
## 6 9.211938 1.086019 7.056770 11.36711 3.375313 0.1891709 2.999909
## b2_upper sigma_hat
## 1 3.292464 6.304427
## 2 3.639503 6.559746
## 3 3.628597 6.609825
## 4 3.377655 6.311905
## 5 3.441463 5.704105
## 6 3.750716 5.662732これで、シミュレーションで得られた数字はsim_1に保存された。
係数の推定値を理解する
予測変数 \(x\) の係数の推定値 b2 の分布を確認してみよう。 私たちは母数\(\beta_2=3\)であることを知っているが、この数値をうまく推定できるだろうか？
hist_b2 <- ggplot(sim_1, aes(x = b2, y = ..density..)) +
geom_histogram(binwidth = .1, color = 'black')
hist_b2 <- hist_b2 +
labs(y = '確率密度', title = 'シミュレーションで得たb2の分布') +
geom_vline(xintercept = 3, color = 'tomato')
print(hist_b2)
このヒストグラムが示すように、線形回帰は平均すると、母数をうまく推定してくれる（不偏性, unbiasedness）。しかし、常に正しい推定をするわけではなく、係数の値を過小推定することもあれば、過大推定することもある。実際の分析では、データセットが1つしかないのが普通であり、自分のデータ分析が係数を「正しく」推定しているとは限らならい。そのために、推定の不確実性を明示することが求められるのである。
標準誤差を理解する
次に、b2の標準誤差 (standar error) をヒストグラムにしてみよう。
hist_se <- ggplot(sim_1, aes(x = b2_se, y = ..density..)) +
geom_histogram(binwidth = .01, color = 'black') +
labs(x = 'b2の標準誤差', y = '確率密度', title = 'b2の標準誤差の分布')
print(hist_se)
このように、標準誤差自体が推定量なので、値はサンプルごとに異なる（分布する）。
標準誤差をseとすると、 \[\frac{\hat{\beta} - \beta}{\mathrm{se}}\] は自由度\(n-k\)（ここでは、100-2=98）の\(t\)分布に従うはずである。 まず、この値をヒストグラムにして、自由度98の\(t\)分布の確率密度曲線を重ね書きしてみよう。
sim_1 <- sim_1 %>% mutate(b2_t = (b2 - beta2) / b2_se)
## t分布の確率密度を求め、true_tという名前のデータフレームに保存する
true_t <- data.frame(x = seq(-4, 4, length = 100)) %>%
mutate(density = dt(x, df = 98))
hist_t <- ggplot(sim_1, aes(x = b2_t, y = ..density..)) +
geom_histogram(binwidth = 0.5, color = 'black') +
geom_line(data = true_t, aes(x = x, y = density), color = 'tomato')
hist_t <- hist_t +
labs(x = expression(frac(hat(beta) - beta, 'se')),
y = '確率密度', title = '標準化されたb2と自由度98のt分布')
print(hist_t)
この図から、\(t\)分布に近い分布であることがわかる。 Q-Qプロットでも確かめてみよう。
qqplot_t <- ggplot(sim_1, aes(sample = b2_t)) +
stat_qq(distribution = qt, dparams = list(df = 98)) +
geom_abline(intercept = 0, slope = 1, color = 'tomato') +
labs(title = 'Q-Q プロット',
x = '自由度n-k=98のt分布',
y = expression(paste('シミュレートした ', frac(hat(beta) - beta, 'se'), ' の分布')))
print(qqplot_t)
Q-Qプロットの点がほぼ一直線に並んでおり、\((\hat{\beta}-\beta)/\mathrm{se}\)が\(t\)分布に従っていることがわかる。ただし、分布の裾では、理論値との乖離が大きいことに注意が必要である。
今回のシミュレーションで得られた標準誤差の平均値は、
## [1] 0.209381である。標準誤差は推定値の標準偏差のはずだが、本当にそうなっているかどうか確認してみよう。
## [1] 0.2144247これで、実際に標準誤差は推定値の標準偏差に（ほぼ）一致することがわかった。
信頼区間を理解する
係数の推定値b2の95パーセント信頼区間を例として考える。 シミュレーションで得られた1つ目の信頼区間は、
## b2_lower b2_upper
## 1 2.388807 3.292464すなわち、[2.39, 3.29] がb2[1] の95パーセント信頼区間である。 この区間は母数である\(\beta_2=3\)を区間内に含んでいる。 したがって、この信頼区間が母数を含む確率は1（100%)である。
同様に、2つ目の信頼区間は、
## b2_lower b2_upper
## 2 2.676187 3.639503であり、これも母数を区間内に含んでいる。
しかし、39番目の信頼区間は、
## b2_lower b2_upper
## 39 3.025416 3.749442であり、これは母数を区間内に含んでいない。 つまり、この信頼区間が母数を含む確率は0である。
シミュレーションで得た1,000個の信頼区間のうち、どの信頼区間が母数を含んでいるか調べてみよう。 母数を信頼区間内に含むのは、信頼区間の下限値が母数以下かつ上限値が母数以上のものである。
この結果、TRUEとなっているものが母数を区間内に捉えているもの、FALSEがそうでないものである。これを表にすると、
## check_ci
## FALSE TRUE
## 49 951と、なる。つまり、1000個の95パーセント信頼区間のうち、951個（95.1%）は母数を区間内に捉え、残りの49個が捉えないということである。このように、一定のデータ生成過程から得られた異なるデートセットに対し、信頼区間を求める作業を繰り返したとき、「得られた信頼区間のうち何パーセントが母数を区間内に含むか」というのが、信頼区間の信頼度である。
これを図にしてみよう。 1000個は多すぎるので、無作為に100個だけ選んで図示する。
# Give an ID to each simulated result
sim_1 <- sim_1 %>% mutate(id = 1:n())
# Put check_ci into the data frame
sim_1$check_ci <- check_ci
# Randomly choose 100 rows of sim values for graph
sim_1_sml <- sim_1 %>% sample_n(100)
# Show simulated ci's graphically
ctplr <- ggplot(sim_1_sml,
aes(x = reorder(id, b2), y = b2,
ymin = b2_lower, ymax = b2_upper,
color = check_ci)) +
geom_hline(yintercept = beta2, linetype = 'dashed') +
geom_pointrange() +
labs(x = 'シミュレーションID', y = 'b2', title = '95％信頼区間') +
scale_color_discrete(name = '母数を', labels = c('含まない', '含む')) +
coord_flip()
print(ctplr)
このように、100個中95個（95パーセント）の95パーセント信頼区間が、母数である3にかかっていることがわかる。