ch12-answers.Rmd

---
title: "Rによる計量政治学：練習問題の解答例"
subtitle: "第12章 回帰分析の前提と妥当性の診断"
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, cache = TRUE)
```

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

```{r load-packs}
library("tidyverse")
if (capabilities("aqua")) { # Macかどうか判定し、Macの場合のみ実行
  theme_set(theme_gray(base_size = 10, base_family = "HiraginoSans-W3"))
}
library("coefplot")
```


## Q12-1

まず、ビールの出荷量に関する分析を実行する（説明は省略する）。
```{r rep-beer-ch11}
#dir.create("data") # dataディレクトリがないなら作る
## データを持っていないならダウンロード
#download.file(url = "https://git.io/fA6Zk",
#              destfile = "data/beer2010.csv")
Beer <- read_csv("data/beer2010.csv")  # データを読み込む
#glimpse(Beer)  # データの中身を確認
## ビールの売り上げを気温に回帰する
fit_beer <- lm(beer ~ temp, data = Beer)
summary(fit_beer)  # 結果を表示
```

この分析の残差プロットを描く。
```{r res-beer}
res_beer <- data_frame(res = fit_beer$residuals,
                       fitted = fit_beer$fitted.values) %>% 
  ggplot(aes(x = fitted, y = res)) +
  geom_point() +
  geom_hline(yintercept = 0) +
  labs(x = "予測値", y = "残差")
print(res_beer)
```

はっきりとしたパタンは見られないものの、隣り合う予測値同士の残差が近いところが多いように見える。また、分散は均一ではないようにも見えるが、一つだけある外れ値を除けば、分散が均一だと言えないこともなさそうだ。標本サイズ（観測数）が12しかないので、はっきりしたことはわからない。

次に、正規QQプロットを描いてみよう。
```{r qq-beer, fig.width = 4, fig.height = 4}
## 残差を標準化する
z_res <- with(fit_beer, (residuals - mean(residuals)) / sd(residuals))
## 図を描く
qq_beer <- data_frame(z_res) %>% 
  ggplot(aes(sample = z_res)) +
    geom_abline(intercept = 0, slope = 1,               # 45度線
                linetype = "dashed", color = "gray") +
    geom_qq() +
    labs(x = "標準正規分布", y = "標準化した残差の分布") +
    xlim(-4, 4) + ylim(-4, 4) +
    coord_fixed()  # 横軸と縦軸を1:1にする
print(qq_beer)
```

点が45度線からずれており、誤差の分布が正規分布に従っていないように見える。標本サイズが小さいので、やはりはっきりしたことはわからないが、このモデルの正しさを積極的に支持することは難しそうである。

## Q12-2

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)
```

2012年衆院選の自民党候補だけを抜き出してデータフレームを作る。
```{r dpj2012}
LDP2012 <- filter(HR, year == 2012, party_jpn == "自民党")
```


### Q11-1-1

得票率 voteshare、選挙費用 expm 、年齢 age の関係を、次の線形モデルで表す。
$$得票率_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$

有意水準を5% (0.05) に設定し、回帰分析でこのモデルの母数である $\beta_0$、$\beta_1$、$\beta_2$ の推定と検定を行う。
```{r fit-q12}
fit_q12 <- lm(voteshare ~ expm + age, data = LDP2012)
summary(fit_q12)
```

包括的な仮説から検定しよう。上で示された結果のうち、最後の行にある$F$統計量の$p値$ を見ると、$0.00783 < 0.05$ だから、有意水準5\%で帰無仮説を棄却する。すなわち、$\beta_1$と$\beta_2$のうち、少なくとも一方は0ではない。言い換えると、（**上で想定したモデルが正しいなら**）選挙費用と年齢のうち少なくとも一方は得票率に影響を与える。

次に、係数を個別に検定しよう。まず、過去の選挙費用 (expm) の係数の$p$値は、$0.000545 < 0.05$ だから、有意水準5\%で帰無仮説を棄却する。すなわち、選挙費用は得票率に影響を与える。$\beta_1$の推定値は約$`r round(coef(fit_q12)[2], 2)`$だから、年齢を一定に保つと、選挙費用が100万円増えるごとに、得票率は平均すると$`r round(coef(fit_q12)[2], 2)`$パーセントポイントずつ上昇すると考えられる。この選挙における自民党候補の選挙費用の標準偏差は、`r round(sd(LDP2012$expm, na.rm = TRUE), 2)`百万円である。選挙費用が400万円（おおよそ1標準偏差分）異なると、$`r round(4 * coef(fit_q12)[2], 2)`$パーセントポイント得票率に差が出るということなので、この効果は実質的にはそれほど重要ではなさそうだ。

年齢 (age) の係数の$p$値は、$0.0139 < 0.05$だから、有意水準5%で帰無仮説を棄却する。すなわち、年齢は得票率に影響を与える。$\beta_2$の推定値は約$`r round(coef(fit_q12)[3], 2)`$だから、過去の当選回数を一定に保つと、年齢が1歳上がるごとに得票率が$`r round(coef(fit_q12)[3], 2)`$パーセントポイントずつ上がることが予測される。
一回り年齢が違うと、$`r round(12 * coef(fit_q12)[3], 2)`$パーセントポイントだけ得票率に差が出るということなので、無視できるほど小さい差ではないが、実質的に大きな効果があると言い切ることは難しい。30歳の候補と60歳の候補を比べると、$`r round(30 * coef(fit_q12)[3], 2)`$パーセントポイントの差になるので、年齢差が大きければ実質的に意味がある違いが出ることがわかる（2012年の自民候補のうち、最年少候補は`r min(LDP2012$age)`歳、最年長候補はは`r max(LDP2012$age)`歳である）。

この分析結果を図にまとめる。
```{r coefplot-q12, fig.width = 8}
coefs_q12 <- coefplot(
  model = fit_q12,    # 図示するモデルを指定
  intercept = FALSE,  # 切片を表示しない
  pointSize = 4,      # 点の大きさを4に
  ## 内側の線を50%信頼区間にする。既定値は1（「点推定値±1標準誤差」）
  innerCI = qt(df = summary(fit_q12)$df[2],
               p = 0.25, lower.tail = FALSE),
  lwdInner = 2.5,     # 内側の線の太さを2.5に
  ## 外側の線を95%信頼区間にする。既定値は2（「点推定値±2標準誤差」）
  outerCI = qt(df = summary(fit_q12)$df[2],
               p = 0.025, lower.tail = FALSE),
  lwdOuter = 0.5,     # 外側の線の太さを0.5に
  ## 説明変数を日本語で表示する
  newNames = list(expm = "選挙費用（100万円）", age = "年齢"),
  xlab = "係数の推定値",
  ylab = "説明変数",
  title = "2012年衆院選における自民党候補の得票率 (%) を予測するモデル"
)
print(coefs_q12)
```

この図が推定結果をまとめている。縦軸にはモデルに含まれる説明変数が並んでいる。横軸は、係数の推定値の大きさを表している。図中の点は、各説明変数の係数の点推定値である。また、水平に引かれた線は、内側の太い線が50パーセント信頼区間、外側の細い線が95パーセント信頼区間である。いずれの変数についても、95パーセント信頼区間がゼロを跨いでいないので、推定された効果は5%の有意水準で統計的に有意であることがわかる。この分析の標本サイズは`r length(summary(fit_q12)$residuals)`、自由度調整済み決定係数は`r round(summary(fit_q12)$adj.r.squared, 2)`である。

### Q12-2-2

上のモデルの残差プロットを作る。
```{r res-ldp-2012}
res_ldp <- data_frame(res = fit_q12$residuals,
                      fitted = fit_q12$fitted.values) %>% 
  ggplot(aes(x = fitted, y = res)) +
  geom_point() +
  geom_hline(yintercept = 0) +
  labs(x = "予測値", y = "残差")
print(res_ldp)
```

特に目立ったパタンはなさそうである。残差の分散は、予測値が小さいとき（予測得票率が40%以下）にそれ以外の場合よりも小さく、誤差の分散は均一ではないかもしれない。

最後に、正規QQプロットで誤差の分布を診断する。
```{r qq-ldp-2012,  fig.width = 4, fig.height = 4}
## 残差を標準化する
z_res <-  with(fit_q12, (residuals - mean(residuals)) / sd(residuals))
## 図を描く
qq_ldp <- data_frame(z_res) %>% 
  ggplot(aes(sample = z_res)) +
    geom_abline(intercept = 0, slope = 1,               # 45度線
                linetype = "dashed", color = "gray") +
    geom_qq(pch = 16, size = 1) +
    labs(x = "標準正規分布", y = "標準化した残差の分布") +
    xlim(-4, 4) + ylim(-4, 4) +
    coord_fixed()  # 横軸と縦軸を1:1にする
print(qq_ldp)
```

分布の端で残差が45度線から外れており、誤差の分布が正規分布ではない可能性が示唆される。