準備

必要なパッケージを読み込む。

pacman::p_load(tidyverse,
               broom)

次に、ggplot2のテーマとフォントの設定を行う。自分好みの設定がある場合は自由に変えてよい。

if (.Platform$OS.type == "windows") { 
  if (require(fontregisterer)) {
    my_font <- "Yu Gothic"
  } else {
    my_font <- "Japan1"
  }
} else if (capabilities("aqua")) {
  my_font <- "HiraginoSans-W3"
} else {
  my_font <- "IPAexGothic"
}

theme_set(theme_gray(base_size = 9,
                     base_family = my_font))

説明のために『Rによる計量政治学』（浅野正彦, 矢内勇生. 2018）で使用されているデータ（hr-data.csv） を使う。 まず、このデータをダウンロードする。

download.file(
  url = "https://raw.githubusercontent.com/yukiyanai/quant-methods-R/master/data_fixed/hr-data.csv",
  dest = "data/hr-data.csv"
)

ダウンロードしたデータを読み込む。

HR <- read_csv("data/hr-data.csv")

小選挙区での当落を表すダミー変数（これを二値の応答変数として扱う）を作る。

HR <- HR %>% 
  mutate(wlsmd = ifelse(wl == "当選", 1, 0))

2009年の結果だけ抜き出し、HR09として保存する。

HR09 <- HR %>% 
  filter(year == 2009)

expが欠測の個体が1つあるので、（とりあえず）除外しておく。

HR09 <- na.omit(HR09)


説明変数が1つのロジスティック回帰

例として、2009年の衆院選データを用い、選挙費用（expm, 単位100万円）で衆院選小選挙区での当落 (wlsmd) を説明するモデルを考える。このモデルを式で表すと、 \[ \mathrm{Pr}(\text{当選}) = \mathrm{logit}^{-1}(\alpha + \beta \text{選挙費用}) \] となる。

まず、2つの変数の関係を図示してみよう。

p <- ggplot(HR09, aes(x = expm, y = wlsmd)) + 
  geom_jitter(size = 1,
              alpha = 1/3,
              width = 0,
              height = 0.05) +
  labs(x = "選挙費用（100万円）",
       y = "小選挙区での当落")
plot(p)

これに（無理やり）通常の回帰直線を当てはめると、

plot(p + geom_smooth(method = "lm"))

となる。

Rでロジスティック回帰分析を行うにはglm() を使う。 この関数は、一般化線形モデル (Generalized Linear Models: GLM) を当てはめるための関数で、ロジスティック回帰以外でも頻繁に使う関数である。 この関数をロジスティック回帰分析に使うには、引数family = binomial(link = "logit") を指定する。

fit1 <- glm(wlsmd ~ expm, 
            data = HR09, 
            family = binomial(link = "logit"))

推定結果を確認する。

tidy(fit1, conf.int = TRUE)
## # A tibble: 2 × 7
##   term        estimate std.error statistic  p.value conf.low conf.high
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
## 1 (Intercept)   -2.01     0.126      -16.0 2.82e-57   -2.26     -1.77 
## 2 expm           0.146    0.0142      10.3 1.14e-24    0.118     0.174

これで、推定はできた。

問題：推定値はどう解釈する？

次に、成功確率の予測値を求めよう。 予測値を求めるには、OLSのときと同様、predict() を使えばよい。 ただし、predict()のデフォルトでは、線形予測子の測定単位で予測値を返す。 ロジスティック回帰の線形予測子\(X\beta\) は対数オッズスケール \((-\infty, \infty)\) なので、 predict(fit1)とすると、成功確率が返ってこない（各自で確かめてみること）。 私たちが欲しいのは確率なので、type = "response" を指定し、応答変数の測定単位で値を得る（ここで、応答 [response] と考えられているのは、当落そのものではなく、当選確率であることに注意）。

HR09 <- HR09 %>% 
  mutate(pred = predict(fit1, type = "response"))

これを図示してみよう。

p_pred <- p +
  geom_line(data = HR09,
            aes(x = expm, y = pred), 
            color = "dodgerblue") +
  labs(x = "選挙費用（100万円）",
       y = "小選挙区での当選確率")
plot(p_pred)

あるいは、ggplot2の機能を利用して、次のようにすることもできる。

p3 <- ggplot(HR09, aes(x = expm, y = wlsmd)) + 
  geom_jitter(size = 1,
              alpha = 1/3,
              width = 0,
              height = 0.05) +
  geom_smooth(method = "glm", 
              method.args = list(family = binomial(link = "logit"))) +
  labs(x = "選挙費用（100万円）",
       y = "小選挙区での当選確率")
print(p3)


説明変数が2つ以上のロジスティック回帰

説明変数が増えても、基本的な分析法は変わらない。 先ほどまでのモデルに、過去の当選回数 (previous) という説明変数を加えてみよう。

ロジスティック回帰を当てはめると、次のような結果が得られる。

fit2 <- glm(wlsmd ~ expm + previous, 
            data = HR09,
            family = binomial(link = "logit"))
tidy(fit2, conf.int = TRUE)
## # A tibble: 3 × 7
##   term        estimate std.error statistic  p.value conf.low conf.high
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
## 1 (Intercept)  -2.03      0.132     -15.4  2.31e-53  -2.30     -1.78  
## 2 expm          0.0500    0.0181      2.77 5.62e- 3   0.0144    0.0853
## 3 previous      0.324     0.0370      8.75 2.06e-18   0.253     0.399

問題：この推定結果をどう解釈する？

この結果を図示したいが、説明変数が複数あるので、先ほどのように簡単に図示することはできない。 選挙費用の効果は統計的に有意でない（ただし、符号は期待通り正）ので、選挙費用を平均値に固定し、過去の当選回数と当選確率の関係を図示してみよう。

選挙費用が平均値のときの当選確率の予測値は次のように求める。

HR09 <- HR09 %>% 
  mutate(pred2 = predict(fit2, 
                         type = "response",
                         list(previous = previous,
                              expm = rep(mean(expm), n()))))

これを図示しよう。

p_pred2 <- ggplot(HR09, aes(x = previous)) + 
  geom_jitter(aes(y = wlsmd), 
              size = 1,
              alpha = 1/3,
              width = 0.1, 
              height = 0.05) +
  geom_line(aes(y = pred2), color = "tomato") +
  labs(x = "当選回数",
       y = "小選挙区での当選確率")
print(p_pred2)

これに信頼区間を加えよう。そのために、まず各予測値の標準誤差を求める。

HR09 <- HR09 %>% 
  mutate(se2 = predict(fit2, type = "response",               
                       list(previous = previous,
                            expm = rep(mean(expm), n())),
                       se.fit = TRUE)$se.fit)

予測値$$2seの範囲を図示する（つまり、ほぼ95%信頼区間を描く）ことにして、その範囲を求める。

HR09 <- HR09 %>% 
  mutate(fit2_lb = pred2 - 2 * se2,
         fit2_ub = pred2 + 2 * se2)

求めた区間を図示しよう。

p_pred2se <- p_pred2 +
  geom_ribbon(data = HR09,
              aes(ymin = fit2_lb, ymax = fit2_ub), 
              alpha = 0.2)
plot(p_pred2se)

ロジスティック回帰の結果に信頼区間を加えることができた。



授業の内容に戻る