比較政治学B

説明変数が二値しかとらないとき（モデル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.843 -12.134  -5.484   8.707  52.016 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  13.9840     0.6221   22.48   <2e-16
## experience   30.8589     0.9832   31.39   <2e-16
## 
## Residual standard error: 16.26 on 1137 degrees of freedom
## Multiple R-squared:  0.4642, Adjusted R-squared:  0.4637 
## F-statistic:   985 on 1 and 1137 DF,  p-value: < 2.2e-16

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

display(fit_1)

## lm(formula = voteshare ~ experience, data = HR09)
##             coef.est coef.se
## (Intercept) 13.98     0.62  
## experience  30.86     0.98  
## ---
## n = 1139, k = 2
## residual sd = 16.26, R-Squared = 0.46

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

この結果を図示しよう。

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

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

Rで議員経験ごとに平均得票率を求め、上の式から求めた予測値と一致するか確かめよう。

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

## # A tibble: 2 x 2
##   experience voteshare
##        <dbl>     <dbl>
## 1          0      14.0
## 2          1      44.8

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

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

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

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

## lm(formula = voteshare ~ expm, data = HR09)
##             coef.est coef.se
## (Intercept) 7.74     0.76   
## expm        3.07     0.10   
## ---
## n = 1124, k = 2
## residual sd = 16.04, R-Squared = 0.48

回帰直線を図示する。

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

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

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

複数の予測変数を使った回帰（モデル3）

次に、得票率を議員経験と選挙費用で説明するモデルを推定する。 lm() で重回帰を行うときは、説明変数を + でつなぐ。

fit_3 <- lm(voteshare ~ experience + expm, data = HR09)
display(fit_3)

## 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 を説明変数として加えてくれるので、これを利用する。

fit_4 <- lm(voteshare ~ experience * expm, data = HR09)
display(fit_4)

## 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

この結果を図示する。

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

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

係数の推定値を理解する

説明変数 \(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\)分布に従っていることがわかる。ただし、分布の裾では、理論値との乖離が大きいことに注意が必要である。

今回のシミュレーションで得られた標準誤差の平均値は、

mean(sim_1$b2_se)

## [1] 0.209381

である。標準誤差は推定値の標準偏差のはずだが、本当にそうなっているかどうか確認してみよう。

sd(sim_1$b2)

## [1] 0.2144247

これで、実際に標準誤差は推定値の標準偏差に（ほぼ）一致することがわかった。

信頼区間を理解する

係数の推定値b2の95パーセント信頼区間を例として考える。 シミュレーションで得られた1つ目の信頼区間は、

sim_1[1, c('b2_lower', 'b2_upper')]

## # A tibble: 1 x 2
##   b2_lower b2_upper
##      <dbl>    <dbl>
## 1     2.39     3.29

すなわち、[2.39, 3.29] がb2[1] の95パーセント信頼区間である。 この区間は母数である\(\beta_2=3\)を区間内に含んでいる。 したがって、この信頼区間が母数を含む確率は1（100%)である。

同様に、2つ目の信頼区間は、

sim_1[2, c('b2_lower', 'b2_upper')]

## # A tibble: 1 x 2
##   b2_lower b2_upper
##      <dbl>    <dbl>
## 1     2.68     3.64

であり、これも母数を区間内に含んでいる。

しかし、39番目の信頼区間は、

sim_1[39, c('b2_lower', 'b2_upper')]

## # A tibble: 1 x 2
##   b2_lower b2_upper
##      <dbl>    <dbl>
## 1     3.03     3.75

であり、これは母数を区間内に含んでいない。 つまり、この信頼区間が母数を含む確率は0である。

シミュレーションで得た1,000個の信頼区間のうち、どの信頼区間が母数を含んでいるか調べてみよう。 母数を信頼区間内に含むのは、信頼区間の下限値が母数以下かつ上限値が母数以上のものである。

check_ci <- sim_1$b2_lower <= beta2 & sim_1$b2_upper >= beta2

この結果、TRUEとなっているものが母数を区間内に捉えているもの、FALSEがそうでないものである。これを表にすると、

table(check_ci)

## 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にかかっていることがわかる。