---
title: "Rによる計量政治学：練習問題の解答例"
subtitle: "第13章 回帰分析の応用"
author: "浅野正彦, 矢内勇生"
date: "2018-11-28"
output:
    html_document:
        theme: united
        highlight: tango
        css: my-markdown.css
        toc: true
        self_contained: true
---

```{r global_option, include = FALSE}
knitr::opts_chunk$set(warning = FALSE, echo = TRUE, message = FALSE,
                      fig.width = 5, fig.height = 4)
```

## パッケージの読み込み

```{r load-packs}
library("tidyverse")
if (capabilities("aqua")) { # Macかどうか判定し、Macの場合のみ実行
  theme_set(theme_gray(base_size = 10, base_family = "HiraginoSans-W3"))
}
library("coefplot")
if(!requireNamespace("knitr")) install.packages("knitr")
library("knitr")
if(!requireNamespace("devtools")) install.packages("devtools")
if(!requireNamespace("patchwork")) devtools::install_github("thomasp85/patchwork")
library("patchwork")
```


## Q13-1

Rds形式の衆院選データ (hr-data.Rds) を読み込む。
手元にない場合はまずダウンロードする。
```{r read-HR-data}
# dir.create("data") # dataディレクトリがない場合は作る
#download.file(url = "https://git.io/fp00p",
#              destfile = "data/hr-data.Rds")
HR <- read_rds("data/hr-data.Rds")
## Rdsファイルの読み込みがうまくいかない場合は以下を実行してCSVファイルを使う
#download.file(url = "https://git.io/fxhQU",
#              destfile = "data/hr-data.csv")
#HR <- read_csv("data/hr-data.csv")
```

正しく読み込めたかどうか確認する。
```{r check-HR}
glimpse(HR)
```

2000年の衆院選だけを抜き出してデータフレームを作る。
```{r dpj2012}
HR2000 <- filter(HR, year == 2000)
```


### Q13-1-1

まず、2000年衆院選の候補者の所属政党 (party\_jpn) を確認する (party を使ってもよい)。
（HTML用の表を作るために、`knitr::kable()` を使う。）
```{r tbl-parties}
# with(HR2000, table(party_jpn))  # R上で確認するだけならこれでよい。
HR2000 %>% 
  with(table(party_jpn)) %>% 
  kable(col.names = c("政党", "候補者数"))
```

無所属とその他を除くと、10の政党（党派）の候補者がいることがわかる。
各政党に所属する候補者の数を見ると、新社会党が1人、無所属の会が9人と人数が少ないので、これらの二つは「その他」にまとめることにする。
```{r parties-group, show = "asis"}
pty <- with(HR2000, ifelse(party_jpn %in% c("新社会党", "無所属の会"), 
                           "その他", party_jpn))
# table(pty)  # R上で確認するならこれでよい
kable(table(pty), col.names = c("政党", "候補者数"))
```

また、回帰分析における参照カテゴリを自民党候補にするため、自民党を第1因子、その他を最後の因子とする因子 (factor) 型の変数を作り、HR2000に加える。
```{r ldp-pos}
(pty_names <- unique(pty))
(pos_ldp <-  which(pty_names == "自民党"))    # 自民党がpty_namesの何番目にあるか確認する
(pos_other <-  which(pty_names == "その他"))  # その他がpty_namesの何番目にあるか確認する
```
自民党は`r pos_ldp`番目、その他は`r pos_other`番目にある。


```{r parties-factor}
HR2000 <- HR2000 %>% 
  mutate(pty_f = factor(pty, levels = c(pty_names[pos_ldp], 
                                        pty_names[-c(pos_ldp, pos_other)],
                                        pty_names[pos_other])))
class(HR2000$pty_f)
#with(HR2000, table(pty_f))  # R上で確認するならこれでよい 
HR2000 %>% 
  with(table(pty_f)) %>% 
  kable(col.names = c("政党", "候補者数"))
```

この変数を使い、以下の線形モデルを仮定して分析する。
$$得票率_i = \beta_0 + \beta_1 選挙費用_i + \sum_{k \in K} \beta_{2k} 政党kのダミー変数_i + \sum_{k \in K} \beta_{3k} 選挙費用_i \cdot 政党kのダミー変数_i + \epsilon_i.$$

ただし、$K=\{$民主党, 公明党, 無所属, 共産党, 自由連合, 保守党, 社民党, 自由党, その他$\}$であり、政党$k$のダミー変数は、候補者$i$の所属政党が$k$のときに1、それ以外のときは0をとる。自民党は参照カテゴリなので、$K$には含まれない。

選挙費用が得票率に与える影響が所属政党によって異なっていたかを確かめるために、上の式で表されるモデルのそれぞれの$k$について、次のような仮説を立てる。

- 帰無仮説：$\beta_{3k} = 0$
- 対立仮説：$\beta_{3k} \neq 0$

有意水準を3% (0.03) に設定し（問題文に指定がないので、自分で決めてよい。慣例では5%だが、常に慣例に従わなければいけないのだろうか）、母数を推定し、この仮説を検定する。得票率は voteshare、選挙費用（100万円）はexpm である。
```{r fit-q13-1-1}
fit_1 <- lm(voteshare ~ expm * pty_f, data = HR2000)
summary(fit_1)
```

$\beta_{3k}$の推定結果に関する情報は、「expm:pty\_f民主党」から「expm:pty\_fその他」までの行に表示されている。
まず、$\beta_{3, 民主}$  の係数の推定値は、「expm:pty\_f民主党」の行の、Estimateの列に表示されている`r round(coef(fit_1)[12], 2)`である。この数値は、ゼロと統計的に意味のある違いを持っているだろうか。これを確かめるために、同じ行の右端の列にある$p$値を確かめる。$p$値は、$0.012217 < 0.03$ なので、$\beta_{3, 民主}$   $= 0$ という帰無仮説を棄却する。よって、選挙費用が得票率に与える影響の大きさは、自民党候補と民主党候補の間で異なるといえる。
自民党候補の場合、選挙費用が得票率に与える影響の大きさは、$`r round(coef(fit_1)[2], 2)`$であるのに対し、民主党候補の場合は、$`r round(coef(fit_1)[2], 2)` + `r round(coef(fit_1)[12], 2)` = `r round(coef(fit_1)[2] + coef(fit_1)[12], 2)`$ である。つまり、自民党候補の場合、他の条件が等しければ、選挙費用を100万円増やすごとに得票率が$`r round(coef(fit_1)[2], 2)`$パーセントポイントずつ上昇することが期待されるが、民主党候補の場合、選挙費用を100万円増やすごとに得票率が`r round(coef(fit_1)[2] + coef(fit_1)[12], 2)`パーセントポイントずつ上昇することが期待される。

他の政党についても、同じように考えればよい。$p$値が0.03（あるいは自分で設定した有意水準）より小さければ、帰無仮説を棄却し、$p$値が0.03以上なら、帰無仮説を棄却しない。例えば、$\beta_{3, 公明}$ の$p$値は$0.31 > 0.03$ なので帰無仮説は棄却されない。すなわち、公明党候補の選挙費用が得票率に与える影響は、自民党候補の選挙費用が得票率に与える影響と異なるという統計的な証拠はない。

$\beta_{3k}$の推定値の数が多いので、キャタピラプロットにまとめて図示する。
```{r catp-beta3, fig.width = 6}
cp1 <- coefplot(
  model = fit_1,      # 図示するモデルを指定
  coefficients = names(fit_1$coefficients[12:20]),  # 指定した係数のみ表示
  sort = "magnitude",                               # 効果が大きい順に並べ替え
  pointSize = 4,      # 点の大きさを4に
  ## 内側の線を50%信頼区間にする。既定値は1（「点推定値±1標準誤差」）
  innerCI = qt(df = summary(fit_1)$df[2],
               p = 0.25, lower.tail = FALSE),
  lwdInner = 2.5,     # 内側の線の太さを2.5に
  ## 外側の線を97%信頼区間にする。既定値は2（「点推定値±2標準誤差」）
  outerCI = qt(df = summary(fit_1)$df[2],
               p = 0.015, lower.tail = FALSE),
  lwdOuter = 0.5,     # 外側の線の太さを0.5に
  xlab = "係数の推定値",
  ylab = "説明変数",
  title = "各党候補者の選挙費用が得票率に与える影響と\n自民党候補の選挙費用が得票率に与える影響の差"
)
print(cp1)
```

この図が$\beta_{3k}$の推定結果をまとめている。縦軸には、どの党について自民党候補との差を見ているかが示されている。
横軸は、$\beta_{3k}$の推定値の大きさを表している。図中の点は、点推定値である。また、水平に引かれた線は、内側の太い線が50パーセント信頼区間、外側の細い線が**97パーセント**信頼区間である（**有意水準を3パーセントに設定した**ので、97パーセントを選んでいる）。公明党と保守党については、97パーセント信頼区間がゼロを跨いでおり、自民党候補との間に統計的に意味のある差が確認できないことがわかる。それ以外の政党については、97パーセント信頼区間がゼロを跨いでいないので、推定された効果は3%の有意水準で統計的に有意であることがわかる。

キャタピラプロットは、選挙費用が得票率に与える影響について「自民党との差」を示しているが、これだけでは各党の候補者の選挙費用が得票率に与えた影響の大きさ自体はわかりにくい。そこで、散布図上に政党別の回帰直線を描く。
```{r ols-parties, fig.width = 8, fig.height = 6}
ols_pty <- ggplot(HR2000, aes(x = expm, y = voteshare, 
                              shape = pty_f, color = pty_f)) +
  geom_smooth(method = "lm", se = FALSE) +
  geom_point() +
  scale_color_discrete(name = "政党") +
  scale_shape_manual(name = "政党", values = 1:10) +
  labs(x = "選挙費用（100万円）", y = "得票率 (%)")
print(ols_pty)
```

得票率と選挙費用の関係が政党ごとに異なることがわかる。しかし、線の数が多すぎて図が見にくいので、各$k \in K$について、自民党と政党$k$を比べる図を作る。

```{r ldp-k}
plt_list <- list()               # 図を保存するためのリスト
pty <- levels(HR2000$pty_f)[-1]  # 自民党の比較対象となる政党
for (i in seq_along(pty)) {      # 比較対象となる政党のループ
  df_plt <- filter(HR2000, pty_f %in% c("自民党", pty[i])) 
  plt <- ggplot(df_plt, aes(x = expm, y = voteshare, 
                            shape = pty_f, color = pty_f)) +
    geom_point() +
    geom_smooth(method = "lm", level = 0.97) +  # 回帰直線と97%信頼区間
    scale_shape_discrete(name = "") +
    scale_color_discrete(name = "") +
    labs(x = "選挙費用", y = "得票率 (%)") +
    theme(legend.position = "top")
  plt_list = c(plt_list, list(plt))
}
```

**patchwork**パッケージを使うと、複数の図を並べて表示することができる（パッケージのインストールと読み込みについては、[パッケージの読み込み](#%e3%83%91%e3%83%83%e3%82%b1%e3%83%bc%e3%82%b8%e3%81%ae%e8%aa%ad%e3%81%bf%e8%be%bc%e3%81%bf) を参照）。


```{r plt-layout, fig.height = 8, fig.width = 8}
plt_list[[1]] + plt_list[[2]] + plt_list[[3]] +
  plt_list[[4]] + plt_list[[5]] + plt_list[[6]] +
  plt_list[[7]] + plt_list[[8]] + plt_list[[9]] +
  plot_layout(ncol = 3, nrow = 3)  # 3行3列に並べて表示する
```

この図は、選挙費用と得票率の関係について、自民党と各政党の候補者を比較している。回帰直線とともに、97パーセント信頼区間が示されている。この図で自民党と比較対象となっている政党それぞれの回帰直線の傾きを見れば、選挙費用が得票率に与える影響の大きさがわかる。例えば、民主党候補と自民党候補の傾きの間には統計的に有意な差があるということは既に述べたが、傾きの差を見ると、それほど大きな違いはないように見える。共産党や社民党の候補は、自民党候補とは大きく異なるようだ。統計的に有意な差が見られなかった公明党や保守党については、自民党候補の傾きとほぼ同じ傾きであることが確認できる。


### Q13-1-2

自民党候補者とそれ以外の候補者を比較するため、自民党候補であることを示すダミー変数 ldp を作る。
```{r ldp-dummy}
HR2000 <- mutate(HR2000, ldp = ifelse(party_jpn == "自民党", 1, 0))
```

この変数を使い、次の線形モデルを仮定して分析する。
$$得票率_i = \beta_0 + \beta_1 選挙費用_i + \beta_{2} 自民党ダミー_i + \beta_{3} 選挙費用_i \times 自民党ダミー_i + \epsilon_i.$$


選挙費用が得票率に与える影響が自民党候補とそれ以外で異なっていたかを確かめるために、次のような仮説を立てる。

- 帰無仮説：$\beta_{3} = 0$
- 対立仮説：$\beta_{3} \neq 0$

有意水準を3% (0.03) に設定し、母数を推定してこの仮説を検定する。
```{r fit-q13-1-2}
fit_2 <- lm(voteshare ~ expm * ldp, data = HR2000)
summary(fit_2)
```
$\beta_{3}$の推定結果に関する情報は、「expm:ldp」の行に表示されている。
まず、$\beta_{3}$  の係数の推定値は、Estimateの列に表示されている$`r round(coef(fit_2)[4], 2)`$である。この数値は、ゼロと統計的に意味のある違いを持っているだろうか。これを確かめるために、同じ行の右端の列にある$p$値を確かめる。$p$値は、2e-16$\approx 0 < 0.03$ なので、$\beta_{3} = 0$ という帰無仮説を棄却する。よって、選挙費用が得票率に与える影響の大きさは、自民党候補とそれ以外の候補の間で異なるといえる。
自民党以外の候補の場合、選挙費用が得票率に与える影響の大きさは、$`r round(coef(fit_2)[2], 2)`$であるのに対し、自民党候補の場合は、$`r round(coef(fit_2)[2], 2)`  `r round(coef(fit_2)[4], 2)` = `r round(coef(fit_2)[2] + coef(fit_2)[4], 2)`$ である。つまり、自民党以外の候補の場合、他の条件が等しければ、選挙費用を100万円増やすごとに得票率が$`r round(coef(fit_2)[2], 2)`$パーセントポイントずつ上昇することが期待されるが、自民党候補の場合、選挙費用を100万円増やすごとに得票率が`r round(coef(fit_2)[2] + coef(fit_2)[4], 2)`パーセントポイントずつ上昇することが期待される。

この結果を図示する。
```{r ols-ldp}
ols_ldp <- ggplot(HR2000, aes(x = expm, y = voteshare, 
                              shape = as.factor(1 - ldp), 
                              color = as.factor(1 - ldp))) +
  geom_smooth(method = "lm", level = 0.97) +    # 回帰直線と97%信頼区間
  scale_shape_discrete(name = "政党", labels = c("自民党", "その他")) +
  scale_color_discrete(name = "政党", labels = c("自民党", "その他")) +
  geom_point() +
  labs(x = "選挙費用（100万円）", y = "得票率 (%)")
print(ols_ldp)
```

この散布図には、自民党候補とそれ以外の候補のそれぞれに、得票率を選挙費用に回帰した直線が引かれている。直線の周りには97パーセント信頼区間が示されている。自民党以外の候補に比べ、自民党候補の回帰直線の傾きが緩やかであることがわかる。つまり、自民党以外の候補と比べると、自民党候補が選挙で費やす費用が得票率に与える影響は小さいことがわかる。


### Q13-1-3

自民党候補だけを抜き出してデータフレームを作る。
```{r df-ldp}
LDP <- filter(HR, party_jpn == "自民党")
## 2000年の自民党候補に限定することも考えられる
# LDP2000 <- filter(HR2000, ldp == 1)
```

このデータを使い、以下の線形モデルを仮定して分析する（交差項を含むモデルを考えたほうがいいが、第14章の内容なので、ここでは割愛する）。
$$
得票率_i = \beta_0 + \beta_1 選挙費用_i + \beta_2 当選回数_i + \epsilon_i.
$$

このモデルの母数について、以下の2組の仮説を立てる。

1. 係数を包括的に検定するための仮説
  - 帰無仮説：$\beta_1 = \beta_2 = 0$
  - 対立仮説：$\beta_1$, $\beta_2$のうち、少なくとも一つは0ではない
2. 係数を個別に検定するための仮説（$k \in \{1, 2\}$）
  - 帰無仮説：$\beta_k = 0$
  - 対立仮説：$\beta_k \neq 0$

回帰分析で、このモデルの母数である $\beta_0$、$\beta_1$、$\beta_2$を推定し、有意水準3パーセントで仮説検定を行う。
```{r fit-3}
fit_3 <- lm(voteshare ~ expm + previous, data = LDP)
summary(fit_3)
```

包括的な仮説から検定しよう。上で示された結果のうち、最後の行にある$F$統計量の$p値$ を見ると、2.2e-16 $= 2.2 \cdot 10^{-16}$未満だから、ほぼ0であり、0.03より小さい。よって、有意水準3\%で帰無仮説を棄却する。すなわち、$\beta_1$と$\beta_2$のうち、少なくとも一つは0ではない。言い換えると、（**上で想定したモデルが正しいなら**）選挙費用と当選回数のうち少なくとも一つは得票率に影響を与える。

次に、係数を個別に検定しよう。まず、選挙費用 (expm) の係数の$p$値は、$0.00541 < 0.03$だから、有意水準3%で帰無仮説を棄却する。すなわち、選挙費用は得票率に影響を与える。
$\beta_1$の推定値は約$`r round(coef(fit_3)[2], 2)`$だから、**当選回数が一定なら**、選挙費用が100万円増えるごとに得票率が$`r round(coef(fit_3)[2], 2)`$パーセントポイントずつ上がることが予測される。自民党候補の選挙費用の標準偏差は約`r round(sd(LDP$expm, na.rm = TRUE))` だが、選挙費用が1標準偏差分増える、すなわち、選挙費用が`r round(sd(LDP$expm, na.rm = TRUE))`百万円増えると、得票率は`r round(coef(fit_3)[2], 2) * round(sd(LDP$expm, na.rm = TRUE))`パーセントポイント上昇する。統計的に有意な結果ではあるが、実質的には有意でない（重要ではない）と判断できる。

次に、過去の当選回数 (previous) の係数の$p$値は、2e-16 $= 2\cdot10^{-16}$未満だから、ほぼ0であり、0.03より小さい。 
よって、有意水準3\%で帰無仮説を棄却する。すなわち、過去の当選回数は得票率に影響を与える。$\beta_2$の推定値は約$`r round(coef(fit_3)[3], 2)`$だから、**選挙費用が一定なら**、過去の当選回数が1回増えるごとに、得票率は平均すると$`r round(coef(fit_3)[3], 2)`$パーセントポイントずつ上昇すると考えられる。過去の当選回数が3回異なると、$`r round(3 *coef(fit_3)[3], 2)`$パーセントポイント得票率に差が出るということなので、この効果は実質的にもある程度意味があると考えられる。

最後に、切片の推定値を確認しよう。$\beta_0$の推定値は約$`r round(coef(fit_3)[1], 2)`$である。これは、選挙費用が0円で、過去の当選回数が0回の候補者の予測得票率が$`r round(coef(fit_3)[1], 2)`$パーセントであることを示している。

### Q13-1-4

[Q13-1-3](#q13-1-3) と同様のモデルを、説明変数を中心化してから推定する。
まず、説明変数を中心化する。
```{r}
LDP <- mutate(LDP,
              c_expm = expm - mean(expm, na.rm = TRUE),
              c_previous = previous - mean(previous))
```

このデータを使い、次の線形モデルを仮定して回帰分析を行う。
$$
得票率_i = \gamma_0 + \gamma_1 中心化した選挙費用_i + \gamma_2 中心化した当選回数_i + \eta_i.
$$

このモデルの母数について、以下の2組の仮説を立てる。

1. 係数を包括的に検定するための仮説
  - 帰無仮説：$\gamma_1 = \gamma_2 = 0$
  - 対立仮説：$\gamma_1$, $\gamma_2$のうち、少なくとも一つは0ではない
2. 係数を個別に検定するための仮説（$k \in \{1, 2\}$）
  - 帰無仮説：$\gamma_k = 0$
  - 対立仮説：$\gamma_k \neq 0$

中心化した説明変数を使って回帰式を推定し、有意水準3%で仮説を検定する。
```{r fit-4}
fit_4 <- lm(voteshare ~ c_expm + c_previous, data = LDP)
summary(fit_4)
```

包括的な仮説から検定しよう。上で示された結果のうち、最後の行にある$F$統計量の$p値$ を見ると、2.2e-16 $= 2.2 \cdot 10^{-16}$未満だから、ほぼ0であり、0.03より小さい。よって、有意水準3\%で帰無仮説を棄却する。すなわち、$\gamma_1$、$\gamma_2$のうち、少なくとも一つは0ではない。言い換えると、（**上で想定したモデルが正しいなら**）選挙費用と当選回数のうち少なくとも一つは得票率に影響を与える。

次に、係数を個別に検定しよう。まず、選挙費用 (expm) の係数の$p$値は、$0.00541 < 0.03$だから、有意水準3%で帰無仮説を棄却する。すなわち、選挙費用は得票率に影響を与える。
$\gamma_1$の推定値は約$`r round(coef(fit_4)[2], 2)`$だから、**当選回数が一定なら**、選挙費用が100万円増えるごとに得票率が$`r round(coef(fit_4)[2], 2)`$パーセントポイントずつ上がることが予測される。自民党候補の選挙費用の標準偏差は`r round(sd(LDP$expm, na.rm = TRUE))` だが、選挙費用が1標準偏差分増える、すなわち、選挙費用が`r round(sd(LDP$expm, na.rm = TRUE))`百万円増えると、得票率は`r round(coef(fit_4)[2], 2) * round(sd(LDP$expm, na.rm = TRUE))`パーセントポイント上昇する。統計的に有意な結果ではあるが、実質的には有意でない（重要ではない）と判断できる。

次に、過去の当選回数 (previous) の係数の$p$値は、2e-16 $= 2\cdot10^{-16}$未満だから、ほぼ0であり、0.03より小さい。 
よって、有意水準3\%で帰無仮説を棄却する。すなわち、過去の当選回数は、得票率に影響を与える。$\gamma_2$の推定値は約$`r round(coef(fit_4)[3], 2)`$だから、**選挙費用が一定なら**、過去の当選回数が1回増えるごとに、得票率は平均すると$`r round(coef(fit_4)[3], 2)`$パーセントポイントずつ上昇すると考えられる。過去の当選回数が3回異なると、$`r round(3 *coef(fit_4)[3], 2)`$パーセントポイント得票率に差が出るということなので、この効果は実質的にもある程度意味があると考えられる。

最後に、切片の推定値を確認しよう。$\gamma_0$の推定値は約$`r round(coef(fit_4)[1], 2)`$である。これは、中心化した選挙費用が0円（すなわち、選挙費用が平均値）で、過去の中心化した過去の当選回数が0回（すなわち、過去の当選回数が平均値）の候補者の予測得票率が$`r round(coef(fit_4)[1], 2)`$パーセントであることを示している。言い換えると、すべての説明変数（中心かする前のもの）が平均値をとるとき、応答変数の予測値は$`r round(coef(fit_4)[1], 2)`$である。


### Q13-1-5

[Q13-1-3](#q13-1-3)の結果と[Q13-1-4](#q13-1-4) の結果を比較すると、実質的な結果はまったく同じであることがわかる。係数の包括的検定の結果も個別的な検定の結果も、説明変数を中心化してもしなくても同じである。また、$\hat{\beta}_1 = \hat{\gamma}_1$、$\hat{\beta}_2 = \hat{\gamma}_2$ である。

しかし、$\hat{\beta}_0 \neq \hat{\gamma}_0$ である。つまり、切片の推定値が異なる。説明変数を中心化しない場合の切片の推定値$\hat{\beta}_0$は、過去の選挙費用が0円で過去の当選回数が0回の候補者の予測得票率を示しているが、そのような候補者はいない。つまり、現実にとってそれほど重要ではない数値が示されている。それに対し、説明変数を中心化した場合の切片の推定値$\hat{\gamma}_0$は、選挙費用と過去の当選回数の両者が平均値の候補者の予測得票率を示している。両方の変数がぴったり平均値に一致する候補者はいないが、この値はデータの中心付近に関する予測値であり、実質的にも意味がある数値である。
