計量経済学応用

説明変数が二値しかとらないとき（モデル1）

得票率（結果変数）を議員経験（説明変数）で説明するモデルを考えよう。 議員経験は、現職または元職の候補者なら1、そうでなければ0をとる二値 (binary) 変数（ダミー変数）である。 このモデルを式で表すと、 \[得票率_i \sim \mathrm{N}(\beta_1 + \beta_2 \cdot 議員経験_i, \sigma^2)\] と、なる。

Rでは、lm() で回帰式を推定する。

fit_1 <- lm(voteshare ~ experience, data = HR09)

これで、fit_1 に推定結果（係数の推定値など）が保存される。

基本的な結果は、summary() で見ることができる。

summary(fit_1)

## 
## Call:
## lm(formula = voteshare ~ experience, data = HR09)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -43.867 -12.072  -5.567   8.583  52.123 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  13.8772     0.6203   22.37   <2e-16
## experience   30.9898     0.9783   31.68   <2e-16
## 
## Residual standard error: 16.19 on 1137 degrees of freedom
## Multiple R-squared:  0.4688, Adjusted R-squared:  0.4684 
## F-statistic:  1004 on 1 and 1137 DF,  p-value: < 2.2e-16

この結果は少し読みにくいので、代わりに brrom:tidy() を利用しよう。

tidy(fit_1)

## # A tibble: 2 x 5
##   term        estimate std.error statistic   p.value
##   <chr>          <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)     13.9     0.620      22.4 3.78e- 92
## 2 experience      31.0     0.978      31.7 2.18e-158

この出力の、estimate の列に係数の推定値 (coefficient estimates) が示されている。 これにより、\(\hat{\beta}_1=\) 13.88, \(\hat{\beta}_2=\) 30.99, \(\hat{\sigma}=\) 16.19 が得られた。 したがって、 \[\widehat{得票率} = 13.88 + 30.99 \cdot 議員経験\] と、なる。

この結果を図示しよう。

p1 <- ggplot(HR09, aes(x = experience, y = voteshare)) +
  scale_x_continuous(breaks = c(0, 1)) +
  geom_jitter(position = position_jitter(width = 0.05), size = 1) +
　geom_smooth(method = 'lm', se = FALSE) + 
  labs(x = '議員経験', y = '得票率（%）')
print(p1 + ggtitle('得票率と議員経験の関係'))

この図に推定の不確実性を示すには、geom_smooth(method = 'lm', se = TRUE) とすればよい（が、 se = TRUE はデフォルトなので、seを指定する必要はない）。 デフォルトでは、95パーセント信頼区間が回帰直線の周りに表示される（この例では、区間が狭すぎてよく見えない）。

p1_ci95 <- p1 + geom_smooth(method = 'lm')
print(p1_ci95 + ggtitle('得票率と議員経験の関係'))

信頼度を変えたいとき、例えば50パーセント信頼区間を表示したいときは、次のようにlevelを指定する。

p1_ci50 <- p1 + geom_smooth(method = 'lm', level = 0.5)
print(p1_ci50 + ggtitle('得票率と議員経験の関係'))

この直線の切片である13.88は、議員経験がない候補者の平均得票率（予測得票率）である。 予測値の式の「議員経験」に0を代入すれば、これは明らかである。 議員経験がある候補者の平均得票率（予測得票率）は、「議員経験」に1を代入することで得られる。 代入してみると、 \(13.88 +30.99 \cdot 1 = 44.87\) となる。

Rで議員経験ごとに平均得票率を求め、上の式から求めた予測値と一致するか確かめよう。 dplyr::group_by() を使うと、指定した変数の値が同じグループを作ることができる。

HR09 %>% 
  group_by(experience) %>%
  summarize(voteshare = mean(voteshare))

## # A tibble: 2 x 2
##   experience voteshare
##        <dbl>     <dbl>
## 1          0      13.9
## 2          1      44.9

このように、予測値は説明変数の値を与えられたときの、結果変数の平均値であることがわかる。

説明変数が連続値をとるとき（モデル2）

同様に、得票率を選挙費用（測定単位：100万円）で説明するモデルは、次のように推定できる。

fit_2 <- lm(voteshare ~ expm, data = HR09)
tidy(fit_2)

## # A tibble: 2 x 5
##   term        estimate std.error statistic   p.value
##   <chr>          <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)     7.74    0.757       10.2 1.61e- 23
## 2 expm            3.07    0.0958      32.1 1.14e-160

回帰直線を図示する。

p2 <- ggplot(HR09, aes(x = expm, y = voteshare)) + 
  geom_point(size = 1) +
  geom_smooth(method = 'lm', se = FALSE) +
  labs(x = '選挙費用（100万円）', y = '得票率（%）') 
print(p2 + ggtitle('得票率と選挙費用の関係'))

95パーセント信頼区間を加える。

p2_ci95 <- p2 + geom_smooth(method = 'lm')
print(p2_ci95 + ggtitle('得票率と選挙費用の関係'))

複数のモデルで回帰分析を実行し、結果を一つの表にまとめたいときは、texreg::screenreg() が便利である（HTMLに出力するなら htmlreg()、LaTeX用には texreg()を使う）。

models <- list(`Model 1` = fit_1,
               `Model 2` = fit_2)
screenreg(models, stars = NULL)

## 
## ===============================
##              Model 1   Model 2 
## -------------------------------
## (Intercept)    13.88      7.74 
##                (0.62)    (0.76)
## experience     30.99           
##                (0.98)          
## expm                      3.07 
##                          (0.10)
## -------------------------------
## R^2             0.47      0.48 
## Adj. R^2        0.47      0.48 
## Num. obs.    1139      1124    
## RMSE           16.19     16.04 
## ===============================

HTMLに出力する場合。

htmlreg(models, stars = NULL, doctype = FALSE)