比較政治学B

4.1.3 世論調査からの予測

選挙結果が保存されている pres08.csv と世論調査結果が保存されている polls08.csv をダウンロードし、プロジェクト内のdataディレクトリに保存する。

download.file(url ="http://qss.princeton.press/student-files/PREDICTION/pres08.csv",
              destfile = "data/pres08.csv")
download.file(url = "http://qss.princeton.press/student-files/PREDICTION/polls08.csv",
              destfile = "data/polls08.csv")

二つのデータセットを読み込んでデータフレームを作り、中身を確認する。

pres08 <- read_csv("data/pres08.csv")
glimpse(pres08)

## Observations: 51
## Variables: 5
## $ state.name <chr> "Alabama", "Alaska", "Arizona", "Arkansas", "Califo...
## $ state      <chr> "AL", "AK", "AZ", "AR", "CA", "CO", "CT", "DC", "DE...
## $ Obama      <int> 39, 38, 45, 39, 61, 54, 61, 92, 62, 51, 47, 72, 36,...
## $ McCain     <int> 60, 59, 54, 59, 37, 45, 38, 7, 37, 48, 52, 27, 62, ...
## $ EV         <int> 9, 3, 10, 6, 55, 9, 7, 3, 3, 27, 15, 4, 4, 21, 11, ...

polls08 <- read_csv("data/polls08.csv")
glimpse(polls08)

## Observations: 1,332
## Variables: 5
## $ state    <chr> "AL", "AL", "AL", "AL", "AL", "AL", "AL", "AL", "AL",...
## $ Pollster <chr> "SurveyUSA-2", "Capital Survey-2", "SurveyUSA-2", "Ca...
## $ Obama    <int> 36, 34, 35, 35, 39, 34, 36, 25, 35, 34, 37, 36, 36, 3...
## $ McCain   <int> 61, 54, 62, 55, 60, 64, 58, 52, 55, 47, 55, 51, 49, 5...
## $ middate  <date> 2008-10-27, 2008-10-15, 2008-10-08, 2008-10-06, 2008...

オバマとマケインの得票差 margin を計算する。

pres08 <- mutate(pres08, margin = Obama - McCain)
polls08 <- mutate(polls08, margin = Obama - McCain)

各州 (+ Washington D.C.) について、最後の世論調査結果を用いてオバマの勝利マージンを予測する。 まず、middate を日付 (Date) クラスに変換し、投票日までの日数を示す DaysToElectionという変数を作る。

polls08 <- polls08 %>% 
  mutate(middate = as.Date(middate),
         DaysToElection = as.Date("2008-11-04") - middate)

次に、州ごとの勝利マージンの予測値を計算する。

st_names <- unique(polls08$state)       # 州の名前を抽出する
poll_pred <- rep(NA, length(st_names))  # 値を記録するベクトル
names(poll_pred) <- st_names
## すべての州でループさせる
for (i in seq_along(st_names)) {
  poll_pred[i] <- polls08 %>%
    filter(state == st_names[i]) %>%  # i番目の州を選ぶ
    filter(DaysToElection == min(DaysToElection)) %>% # 最後の調査
    with(mean(margin)) 
}

最後の世論調査の誤差（poll_predに州の名前が付いているので、errorsにも自動的に州の名前が付く）。

errors <- pres08$margin - poll_pred

平均予測誤差。

mean(errors)

## [1] 1.062092

二乗平均平方根誤差 (root mean squared error: RMSE)。

sqrt(mean(errors^2))

## [1] 5.90894

予測誤差のヒストグラム。

hist_errors <- data_frame(errors = errors) %>%
  ggplot(aes(x = errors, y = ..density..)) +
    geom_histogram(binwidth = 5, boundary = 0, color = "black") +
    geom_vline(xintercept = mean(errors), 
               linetype = "dashed", color = "red") +
    labs(x = "オバマの予測マージンにおける誤差\n（パーセントポイント）",
         y = "確率密度", title = "世論調査からの予測誤差") 
print(hist_errors)

世論調査結果を横軸に、実際の選挙結果を縦軸にしたプロットを作る。

pres08$poll_pred <- poll_pred  # 世論調査による予測値を選挙結果データに加える
real_v_pred <- ggplot(pres08, aes(x = poll_pred, y = margin)) +
  xlim(-50, 100) + ylim(-50, 100) +  # 作図範囲の指定
  geom_vline(xintercept = 0, color = "gray") +  # x=0の縦線
  geom_hline(yintercept = 0, color = "gray") +  # y=0の横線
  geom_abline(intercept = 0, slope = 1, linetype = "dashed") + # 45度線
  geom_text(aes(label = state)) + # stateの値をラベルとしてプロット
  coord_fixed() +  # 縦横比を揃える
  labs(x = "世論調査結果", y = "実際の選挙結果")
print(real_v_pred)

どの州の世論調査が間違えたのか？また、それらの州の実際のマージンはいくつだったか？

pres08 %>%
  filter(sign(poll_pred) != sign(margin)) %>%
  select(state, margin)

## # A tibble: 3 x 2
##   state margin
##   <chr>  <int>
## 1 IN         1
## 2 MO        -1
## 3 NC         1

実際にオバマがえた選挙人票の総数。

pres08 %>%
  filter(margin > 0) %>%
  with(sum(EV))

## [1] 364

世論調査に基づく予測。

pres08 %>%
  filter(poll_pred > 0) %>%
  with(sum(EV))

## [1] 349

選挙期間を通じた世論調査の結果を保存した pollsUS08.csv をダウンロードし、dataディレクトリに保存する。

download.file(url = "http://qss.princeton.press/student-files/PREDICTION/pollsUS08.csv",
              destfile = "data/pollsUS08.csv")

このデータセットを読み込み、中身を確認する。

pollsUS08 <- read_csv("data/pollsUS08.csv")
glimpse(pollsUS08)

## Observations: 526
## Variables: 4
## $ Pollster <chr> "IBD/TIPP", "GWU (Lake/Tarrance)", "Diageo/Hotline", ...
## $ McCain   <int> 48, 51, 43, 44, 44, 45, 49, 42, 44, 48, 44, 43, 40, 4...
## $ Obama    <int> 36, 39, 38, 46, 47, 47, 42, 48, 44, 48, 41, 43, 44, 4...
## $ middate  <date> 2007-07-04, 2007-07-11, 2007-07-14, 2007-07-19, 2007...

投票日までの日数を計算する。

pollsUS08 <- pollsUS08 %>%
  mutate(middate = as.Date(middate),
         DaysToElection = as.Date("2008-11-04") - middate)

予測を保存するための空のベクトルを二つ用意する。

Obama_pred <- McCain_pred <- rep(NA, 90)

ループを使って、7日分の世論調査の平均値を、対象日を動かして計算する。

for (i in 1:90) {
  week_data <- pollsUS08 %>% 
    filter(DaysToElection <= (90 - i + 7),
           DaysToElection > (90 - i))
  Obama_pred[i] <- mean(week_data$Obama)
  McCain_pred[i] <- mean(week_data$McCain)  
}

投票日90日前から前日までの各候補者に対する支持の変化をプロットする。

ts_support <- data_frame(O = Obama_pred, M = McCain_pred,
                         DaysToElection = 90:1) %>%
  ggplot(aes(x = DaysToElection)) +
  geom_point(aes(y = O), color = "blue", pch = 1) +
  geom_point(aes(y = M), color = "red", pch = 1) +
  geom_vline(xintercept = 0, color = "gray") +
  geom_point(aes(x = 0, y = 52.93), pch = 19, color = "blue") +
  geom_point(aes(x = 0, y = 45.65), pch = 19, color = "red") +
  geom_text(aes(x = 80, y = 48, label = "Obama"), color = "blue") +
  geom_text(aes(x = 80, y = 41, label = "McCain"), color = "red") +
  xlim(0, 90) + ylim(40, 60) +
  scale_x_reverse(breaks = seq(0, 90, by = 10)) +
  labs(x = "選挙までの日数",
       y = "候補者への支持\n（パーセントポイント）")
print(ts_support)

4.2.1 顔の見た目と選挙結果

顔の見た目による能力評価と選挙結果に関するデータを保存した face.csv をダウンロードし、dataディレクトリに保存する。

download.file(url = "http://qss.princeton.press/student-files/PREDICTION/face.csv",
              destfile = "data/face.csv")

ダウンロードしたデータセットを読み込み、中身を確認する。

face <- read_csv("data/face.csv")
glimpse(face)

## Observations: 119
## Variables: 10
## $ year    <int> 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, ...
## $ state   <chr> "CA", "DE", "FL", "GA", "HI", "IN", "MA", "MD", "ME", ...
## $ winner  <chr> "Feinstein", "Carper", "Nelson", "Miller", "Akaka", "L...
## $ loser   <chr> "Campbell", "Roth", "McCollum", "Mattingly", "Carroll"...
## $ w.party <chr> "D", "D", "D", "D", "D", "R", "D", "D", "R", "D", "D",...
## $ l.party <chr> "R", "R", "R", "R", "R", "D", "R", "R", "D", "R", "R",...
## $ d.comp  <dbl> 0.5645676, 0.3419122, 0.6123680, 0.5415328, 0.6802323,...
## $ r.comp  <dbl> 0.4354324, 0.6580878, 0.3876320, 0.4584672, 0.3197677,...
## $ d.votes <int> 5790154, 181387, 2987644, 1390428, 251130, 684242, 187...
## $ r.votes <int> 3779325, 142683, 2703608, 933698, 84657, 1419629, 3347...

民主党の得票率 d_share、共和党の得票率 r_share と民主党の勝利マージン diff_share という三つの変数を作る。

face <- face %>%
  mutate(d_share = d.votes / (d.votes + r.votes),
         r_share = r.votes / (d.votes + r.votes),
         diff_share = d_share - r_share)

民主党の勝利マージンと民主党候補に関する見た目の能力の散布図を作る。

scat_face <- ggplot(face, aes(x = d.comp, y = diff_share)) +
  geom_point(color = ifelse(face$w.party == "R", "red", "blue")) +
  xlim(0, 1) + ylim(-1, 1) +
  labs(x = "民主党候補の見た目の能力スコア", 
       y = "民主党候補の得票率マージン",
       title = "見た目の能力と得票率") +
  coord_fixed(ratio = .5)
print(scat_face)

4.2.2 相関と散布図

見た目の能力と勝利マージン（実際の選挙での強さ）の相関係数を求める。

with(face, cor(d.comp, diff_share))

## [1] 0.4327743

4.2.3 最小二乗法

勝利マージンを見た目の能力スコアに回帰する。言い換えると、見た目の能力スコアを使って勝利マージンを予測するモデルを推定する。

fit <- lm(diff_share ~ d.comp, data = face)

推定結果を確認する。

summary(fit)

## 
## Call:
## lm(formula = diff_share ~ d.comp, data = face)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.67487 -0.16600  0.01399  0.17741  0.74297 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.31223    0.06596  -4.733 6.24e-06
## d.comp       0.66038    0.12718   5.193 8.85e-07
## 
## Residual standard error: 0.2665 on 117 degrees of freedom
## Multiple R-squared:  0.1873, Adjusted R-squared:  0.1803 
## F-statistic: 26.96 on 1 and 117 DF,  p-value: 8.854e-07

この結果は、勝利マージンと見た目のスコアの間に以下の関係がある（かもしれない）ことを示す。 \[\widehat{\mathrm{勝利マージン}} = -0.312 + 0.660 \cdot \mathrm{見た目の能力スコア}. \]

回帰分析の推定結果に出てくる \(\hat{}\) は「ハット」と読み（\(\hat{y}\)は\(y\)ハット）、推定値または予測値であることを示す。

回帰分析で予測される変数のことを、結果変数 (outcome varible) または応答変数 (response variable) と呼ぶ。上の分析では、民主党候補の勝利マージンが結果変数である。 それに対し、予測に用いる変数のことを説明変数 (explanatory variable[s]) と呼ぶ。 上の分析では、見た目の勝利スコアが説明変数である。回帰分析では、結果変数を説明変数に回帰する（逆ではないので注意）。回帰分析のうち、説明変数が一つのものを単回帰 (simple regression)、説明変数が二つ以上のものを重回帰 (multiple regression) と呼ぶ。上の分析は単回帰である。

単回帰の場合、ggplot2 を使うと回帰直線を散布図に簡単に上書きできる。具体的には、geom_smooth(method = "lm") を使えばよい。 散布図は先ほど作ったので、先ほどの散布図に回帰直線を加える。

scat_fit <- scat_face +
  geom_smooth(method = "lm", se = FALSE, color = "black")
print(scat_fit)

4.2.4 平均への回帰

省略

4.2.5 Rにおけるデータの結合

2012年の合衆国大統領戦のデータファイル pres12.csv をダウンロードし、dataディレクトリに保存する。

download.file(url = "http://qss.princeton.press/student-files/PREDICTION/pres12.csv",
              destfile = "data/pres12.csv")

データを読み込み、中身を確認する。

pres12 <- read_csv("data/pres12.csv")
glimpse(pres12)

## Observations: 51
## Variables: 4
## $ state  <chr> "AL", "AK", "AZ", "AR", "CA", "CO", "CT", "DE", "DC", "...
## $ Obama  <int> 38, 41, 45, 37, 60, 51, 58, 59, 91, 50, 45, 71, 33, 58,...
## $ Romney <int> 61, 55, 54, 61, 37, 46, 41, 40, 7, 49, 53, 28, 65, 41, ...
## $ EV     <int> 9, 3, 11, 6, 55, 9, 7, 3, 3, 29, 16, 4, 4, 20, 11, 6, 6...

pres08もざっと見る。

glimpse(pres08)

## Observations: 51
## Variables: 7
## $ state.name <chr> "Alabama", "Alaska", "Arizona", "Arkansas", "Califo...
## $ state      <chr> "AL", "AK", "AZ", "AR", "CA", "CO", "CT", "DC", "DE...
## $ Obama      <int> 39, 38, 45, 39, 61, 54, 61, 92, 62, 51, 47, 72, 36,...
## $ McCain     <int> 60, 59, 54, 59, 37, 45, 38, 7, 37, 48, 52, 27, 62, ...
## $ EV         <int> 9, 3, 10, 6, 55, 9, 7, 3, 3, 27, 15, 4, 4, 21, 11, ...
## $ margin     <int> -21, -21, -9, -20, 24, 9, 23, 85, 25, 3, -5, 45, -2...
## $ poll_pred  <dbl> -25.0, -19.0, -2.5, -7.0, 24.0, 7.0, 25.0, 69.0, 30...

二つのデータフレームを結合する。教科書の例と同じように、列と列を結合する（言い換えると、データフレームを横に並べてくっつける）。列の結合には、dplyrパッケージの xxx_join()を使う。xxxの部分は目的によって異なり、二つのデータフレームの観察対象（行）の内容異なる場合の対処法に応じて、以下のような関数がある（この他に、semi_join()もあるが割愛する。）。

• left_join()：左側のデータフレームに存在する観察のみ残す
• right_join()：右側のデータフレームに存在する観察のみ残す
• full_join()：左右どちらかのデータフレームに存在する観察はすべて残す
• inner_join()：両方のデータフレームに存在する観察のみ残す

したがって、一般的には inner_join()の結果がもっとも小さく、full_join() の結果がもっとも大きい。二つのデータフレームに存在する観察対象がまったく同じ場合には、どれを使っても同じ結果が出る。二つのデータフレームが同じかどうかを判定する変数は、引数 by で指定する。

pres08 とpres12を結合してみよう。

pres <- full_join(pres08, pres12, by = "state")

結合されたデータフレームの中身を確認する。

glimpse(pres)

## Observations: 51
## Variables: 10
## $ state.name <chr> "Alabama", "Alaska", "Arizona", "Arkansas", "Califo...
## $ state      <chr> "AL", "AK", "AZ", "AR", "CA", "CO", "CT", "DC", "DE...
## $ Obama.x    <int> 39, 38, 45, 39, 61, 54, 61, 92, 62, 51, 47, 72, 36,...
## $ McCain     <int> 60, 59, 54, 59, 37, 45, 38, 7, 37, 48, 52, 27, 62, ...
## $ EV.x       <int> 9, 3, 10, 6, 55, 9, 7, 3, 3, 27, 15, 4, 4, 21, 11, ...
## $ margin     <int> -21, -21, -9, -20, 24, 9, 23, 85, 25, 3, -5, 45, -2...
## $ poll_pred  <dbl> -25.0, -19.0, -2.5, -7.0, 24.0, 7.0, 25.0, 69.0, 30...
## $ Obama.y    <int> 38, 41, 45, 37, 60, 51, 58, 91, 59, 50, 45, 71, 33,...
## $ Romney     <int> 61, 55, 54, 61, 37, 46, 41, 7, 40, 49, 53, 28, 65, ...
## $ EV.y       <int> 9, 3, 11, 6, 55, 9, 7, 3, 3, 29, 16, 4, 4, 20, 11, ...

左右のデータフレームで異なる名前の変数を利用して結合する場合は以下のようにする。

names(pres12)[1] <- "state_abb"  # 練習のため、変数名を変える
pres <- full_join(pres08, pres12, by = c("state" = "state_abb"))

結合されたデータフレームの中身を確認する。

glimpse(pres)

## Observations: 51
## Variables: 10
## $ state.name <chr> "Alabama", "Alaska", "Arizona", "Arkansas", "Califo...
## $ state      <chr> "AL", "AK", "AZ", "AR", "CA", "CO", "CT", "DC", "DE...
## $ Obama.x    <int> 39, 38, 45, 39, 61, 54, 61, 92, 62, 51, 47, 72, 36,...
## $ McCain     <int> 60, 59, 54, 59, 37, 45, 38, 7, 37, 48, 52, 27, 62, ...
## $ EV.x       <int> 9, 3, 10, 6, 55, 9, 7, 3, 3, 27, 15, 4, 4, 21, 11, ...
## $ margin     <int> -21, -21, -9, -20, 24, 9, 23, 85, 25, 3, -5, 45, -2...
## $ poll_pred  <dbl> -25.0, -19.0, -2.5, -7.0, 24.0, 7.0, 25.0, 69.0, 30...
## $ Obama.y    <int> 38, 41, 45, 37, 60, 51, 58, 91, 59, 50, 45, 71, 33,...
## $ Romney     <int> 61, 55, 54, 61, 37, 46, 41, 7, 40, 49, 53, 28, 65, ...
## $ EV.y       <int> 9, 3, 11, 6, 55, 9, 7, 3, 3, 29, 16, 4, 4, 20, 11, ...

右側のstate_abb は左側のstate とマッチされ、残っていないことがわかる。

オバマの得票率の\(z\)得点（標準化された値）を求める。

pres <- pres %>%
  mutate(Obama2008_z = scale(Obama.x),
         Obama2012_z = scale(Obama.y))

オバマの2012年の標準化得票率を2008年の標準化得票率に回帰する。 ｀

fit1 <- lm(Obama2012_z ~ Obama2008_z, data = pres)
summary(fit1)

## 
## Call:
## lm(formula = Obama2012_z ~ Obama2008_z, data = pres)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.49061 -0.10105 -0.00235  0.12036  0.50853 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.470e-17  2.563e-02    0.00        1
## Obama2008_z  9.834e-01  2.589e-02   37.99   <2e-16
## 
## Residual standard error: 0.1831 on 49 degrees of freedom
## Multiple R-squared:  0.9672, Adjusted R-squared:  0.9665 
## F-statistic:  1443 on 1 and 49 DF,  p-value: < 2.2e-16

散布図と回帰直線を示す。

plt_Obama <- ggplot(pres, aes(x = Obama2008_z, y = Obama2012_z)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  xlim(-4, 4) + ylim(-4, 4) +
  labs(x = "オバマの標準化得票率（2008年）",
       y = "オバマの標準化得票率（2008年）")
print(plt_Obama)

オバマの2008年の得票率が下位25%だった州のうち、2012年の選挙の得票率のほうが高かった州の割合を計算する。

pres %>%
  filter(Obama2008_z <= quantile(Obama2008_z, prob = 0.25)) %>%
  with(mean(Obama2012_z > Obama2008_z))

## [1] 0.5714286

オバマの2008年の得票率が上位25%だった州のうち、2012年の選挙の得票率のほうが高かった州の割合を計算する。

pres %>%
  filter(Obama2008_z > quantile(Obama2008_z, prob = 0.75)) %>%
  with(mean(Obama2012_z > Obama2008_z))

## [1] 0.4615385

4.2.6 モデルの当てはまり

フロリダ州における大統領選挙結果を保存したデータファイル florida.csv をダウンロードし、dataディレクトリに保存する。

download.file(url = "http://qss.princeton.press/student-files/PREDICTION/florida.csv", 
              destfile = "data/florida.csv")

データを読み込み中身を確認する。

florida <- read_csv("data/florida.csv")
glimpse(florida)

## Observations: 67
## Variables: 7
## $ county     <chr> "Alachua", "Baker", "Bay", "Bradford", "Brevard", "...
## $ Clinton96  <int> 40144, 2273, 17020, 3356, 80416, 320736, 1794, 2712...
## $ Dole96     <int> 25303, 3684, 28290, 4038, 87980, 142834, 1717, 2783...
## $ Perot96    <int> 8072, 667, 5922, 819, 25249, 38964, 630, 7783, 7244...
## $ Bush00     <int> 34124, 5610, 38637, 5414, 115185, 177323, 2873, 354...
## $ Gore00     <int> 47365, 2392, 18850, 3075, 97318, 386561, 2155, 2964...
## $ Buchanan00 <int> 263, 73, 248, 65, 570, 788, 90, 182, 270, 186, 122,...

ブキャナンの2000年の得票数を1996年のペローの得票数に回帰する。

fit2 <- lm(Buchanan00 ~ Perot96, data = florida)
summary(fit2)

## 
## Call:
## lm(formula = Buchanan00 ~ Perot96, data = florida)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -612.74  -65.96    1.94   32.88 2301.66 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  1.34575   49.75931   0.027    0.979
## Perot96      0.03592    0.00434   8.275 9.47e-12
## 
## Residual standard error: 316.4 on 65 degrees of freedom
## Multiple R-squared:  0.513,  Adjusted R-squared:  0.5055 
## F-statistic: 68.48 on 1 and 65 DF,  p-value: 9.474e-12

決定係数 (\(R^2\)) は Multiple R-squared として表示されている。 重回帰の場合は、決定係数の代わりに自由度調整済み決定係数 (Adjusted \(R^2\)) を使う。 決定係数は重要ではないので、これ以上は説明しない（教科書を参照）。

残差プロットを作成する。

# データフレームに予測値と残差を加える
florida <- florida %>% 
  mutate(fitted2 = fitted(fit2),
         resid2 = resid(fit2))
# 作図
res_plt2 <- florida %>%
  ggplot(aes(x = fitted2, y = resid2)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "gray") +
  labs(x = "予測値", y = "残差")
print(res_plt2)

外れ値はどの郡 (county) か、すなわち残差が一番大きい郡を調べる。which.max() と使うと、与えた変数の値が一番大きくなる位置を教えてくれるので、これを使うと簡単である。

florida$county[which.max(florida$resid2)]

## [1] "PalmBeach"

Palm Beach 郡を除いた分析を行う。新しいデータフレームを作らなくても、filter()を実行してPalm Beachを除外した後に、パイプでデータフレームを渡せばよい。第1引数以外に渡すので、“.”（ドット）で渡す位置を明示する。

fit3 <- florida %>%
  filter(county != "PalmBeach") %>%
  lm(Buchanan00 ~ Perot96, data = .)  # .の位置にデータフレームを渡す
summary(fit3)

## 
## Call:
## lm(formula = Buchanan00 ~ Perot96, data = .)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -206.70  -43.51  -16.02   26.92  269.03 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) 45.841933  13.892746    3.30  0.00158
## Perot96      0.024352   0.001273   19.13  < 2e-16
## 
## Residual standard error: 87.75 on 64 degrees of freedom
## Multiple R-squared:  0.8512, Adjusted R-squared:  0.8488 
## F-statistic:   366 on 1 and 64 DF,  p-value: < 2.2e-16

残差プロットを作る。 残差プロットを作成する。

res_plt3 <- florida %>% 
  filter(county != "PalmBeach") %>% 
  mutate(fitted3 = fitted(fit3),
         resid3 = residuals(fit3)) %>% 
  ggplot(aes(x = fitted3, y = resid3)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "gray") +
  labs(x = "予測値", y = "残差")
print(res_plt3)

散布図と二つの回帰直線を描く。

scat_fit2_3 <- ggplot(florida, aes(x = Perot96, y = Buchanan00)) +
  geom_smooth(method = "lm", se = FALSE, 
              color = "black", linetype = "dashed") +
  geom_abline(intercept = coef(fit3)[1], slope = coef(fit3)[2]) +
  geom_point() +
  geom_text(aes(x = 30000, y = 3200, label = "Palm Beach")) +
  geom_text(aes(x = 20000, y = 1250, label = "Palm Beachを含んだ回帰直線"),
            family = "HiraginoSans-W3") +
  geom_text(aes(x = 32000, y = 500, label = "Palm Beachを除いた\n回帰直線"),
            family = "HiraginoSans-W3") +
  labs(x = "1996年のペローの得票数",
       y = "2000年のブキャナンの得票数")
print(scat_fit2_3)

4.3.1 ランダム化実験

女性議長のランダム割当実験のデータファイル women.csv をダウンロードし、dataディレクトリに保存する。

download.file(url = "http://qss.princeton.press/student-files/PREDICTION/women.csv",
              destfile = "data/women.csv")

データセットを読み込み、中身を確認する。

women <- read_csv("data/women.csv")
glimpse(women)

## Observations: 322
## Variables: 6
## $ GP         <int> 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, ...
## $ village    <int> 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, ...
## $ reserved   <int> 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, ...
## $ female     <int> 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, ...
## $ irrigation <int> 0, 5, 2, 4, 0, 0, 4, 0, 0, 0, 0, 4, 0, 0, 4, 5, 0, ...
## $ water      <int> 10, 0, 2, 31, 0, 0, 7, 12, 28, 0, 23, 12, 0, 0, 41,...

割り当てのあるGPと割り当てのないGPで女性議長の割合を計算する。

women %>% 
  group_by(reserved) %>%
  summarize(ratio = mean(female))

## # A tibble: 2 x 2
##   reserved  ratio
##      <int>  <dbl>
## 1        0 0.0748
## 2        1 1

飲用水設備に女性議長が与える因果効果

women %>% 
  group_by(reserved) %>%
  summarize(water = mean(water)) %>%
  with(water[2] - water[1])  

## [1] 9.252423

灌漑設備に女性議長が与える因果効果

women %>% 
  group_by(reserved) %>%
  summarize(irrigation = mean(irrigation)) %>%
  with(irrigation[2] - irrigation[1])  

## [1] -0.3693319

回帰分析で平均因果効果を推定する。

fit_water <- lm(water ~ reserved, data = women)
summary(fit_water)

## 
## Call:
## lm(formula = water ~ reserved, data = women)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -23.991 -14.738  -7.865   2.262 316.009 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   14.738      2.286   6.446 4.22e-10
## reserved       9.252      3.948   2.344   0.0197
## 
## Residual standard error: 33.45 on 320 degrees of freedom
## Multiple R-squared:  0.01688,    Adjusted R-squared:  0.0138 
## F-statistic: 5.493 on 1 and 320 DF,  p-value: 0.0197

fit_irrigation <- lm(irrigation ~ reserved, data = women)
summary(fit_irrigation)

## 
## Call:
## lm(formula = irrigation ~ reserved, data = women)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.388 -3.388 -3.019 -1.019 86.612 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   3.3879     0.6498   5.214 3.33e-07
## reserved     -0.3693     1.1220  -0.329    0.742
## 
## Residual standard error: 9.506 on 320 degrees of freedom
## Multiple R-squared:  0.0003385,  Adjusted R-squared:  -0.002785 
## F-statistic: 0.1084 on 1 and 320 DF,  p-value: 0.7422

4.3.2 重回帰モデル

社会的プレッシャーと投票率に関するデータセット social.csv （第2章で使ったので、ダウロード済みのはず）を読み込み、中身を確認する。

social <- read_csv("data/social.csv")
glimpse(social)

## Observations: 305,866
## Variables: 6
## $ sex         <chr> "male", "female", "male", "female", "female", "mal...
## $ yearofbirth <int> 1941, 1947, 1951, 1950, 1982, 1981, 1959, 1956, 19...
## $ primary2004 <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 0,...
## $ messages    <chr> "Civic Duty", "Civic Duty", "Hawthorne", "Hawthorn...
## $ primary2006 <int> 0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 1, 1, 1,...
## $ hhsize      <int> 2, 2, 3, 3, 3, 3, 3, 3, 2, 2, 1, 2, 2, 1, 2, 2, 1,...

処置の違いは、messages 変数に記録されている。

with(social, table(messages))

## messages
## Civic Duty    Control  Hawthorne  Neighbors 
##      38218     191243      38204      38201

比較の基準となるグループが Control になるように、messages 変数を変換する。

social <- social %>% 
  mutate(messages = factor(messages, 
                           levels = c("Control", "Civic Duty",
                                      "Hawthorne", "Neighbors")))
levels(social$messages)

## [1] "Control"    "Civic Duty" "Hawthorne"  "Neighbors"

指定した順番通りの因子型変数ができた。

この因子型変数に2006年の予備選挙での投票を回帰する。

fit_vt <- lm(primary2006 ~ messages, data = social)
summary(fit_vt)

## 
## Call:
## lm(formula = primary2006 ~ messages, data = social)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.3780 -0.2966 -0.2966  0.6776  0.7034 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)
## (Intercept)        0.296638   0.001058 280.393  < 2e-16
## messagesCivic Duty 0.017899   0.002592   6.905 5.03e-12
## messagesHawthorne  0.025736   0.002593   9.927  < 2e-16
## messagesNeighbors  0.081310   0.002593  31.360  < 2e-16
## 
## Residual standard error: 0.4627 on 305862 degrees of freedom
## Multiple R-squared:  0.003283,   Adjusted R-squared:  0.003273 
## F-statistic: 335.8 on 3 and 305862 DF,  p-value: < 2.2e-16

基準カテゴリが教科書と異なるので、教科書とは異なる数字が表示される。しかし、実質的には同じ結果である。例えば、教科書ではCivic Dutyが基準カテゴリでControlグループの係数が\(-0.017899\) だが、ここではControlが基準カテゴリで Civic Dutyグループの係数が\(0.017899\)だから、Civic DutyとControlの差が\(0.017899\)で同じである。

同様に、指標変数（特定のグループに所属していることを示すダミー変数）を用いて分析しなおしてみよう。ただし、ここでも基準カテゴリはControlにする。

social <- social %>% 
  mutate(CD = ifelse(messages == "Civic Duty", 1, 0),
         Hawthorne = ifelse(messages == "Hawthorne", 1, 0),
         Neighbors = ifelse(messages == "Neighbors", 1, 0))
fit_vt2 <- lm(primary2006 ~ CD + Hawthorne + Neighbors, data = social)
summary(fit_vt2)

## 
## Call:
## lm(formula = primary2006 ~ CD + Hawthorne + Neighbors, data = social)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.3780 -0.2966 -0.2966  0.6776  0.7034 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.296638   0.001058 280.393  < 2e-16
## CD          0.017899   0.002592   6.905 5.03e-12
## Hawthorne   0.025736   0.002593   9.927  < 2e-16
## Neighbors   0.081310   0.002593  31.360  < 2e-16
## 
## Residual standard error: 0.4627 on 305862 degrees of freedom
## Multiple R-squared:  0.003283,   Adjusted R-squared:  0.003273 
## F-statistic: 335.8 on 3 and 305862 DF,  p-value: < 2.2e-16

先ほどと同じ結果が得られた。

予測値を計算するため、重複しない “messages” の値からなるデータフレームを作る。

unique_messages <- data_frame(messages = unique(social$messages))
unique_messages

## # A tibble: 4 x 1
##   messages  
##   <fct>     
## 1 Civic Duty
## 2 Hawthorne 
## 3 Control   
## 4 Neighbors

このデータフレームから観察ごとの予測を行う。

predict(fit_vt, newdata = unique_messages)

##         1         2         3         4 
## 0.3145377 0.3223746 0.2966383 0.3779482

messagesの値ごとにprimary2006の標本平均を求める。

social %>% 
  group_by(messages) %>%
  summarize(mean_turnout = mean(primary2006))

## # A tibble: 4 x 2
##   messages   mean_turnout
##   <fct>             <dbl>
## 1 Control           0.297
## 2 Civic Duty        0.315
## 3 Hawthorne         0.322
## 4 Neighbors         0.378

切片なしの回帰モデル。

fit_noint <- lm(primary2006 ~ -1 + messages, data = social)
summary(fit_noint)

## 
## Call:
## lm(formula = primary2006 ~ -1 + messages, data = social)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.3780 -0.2966 -0.2966  0.6776  0.7034 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)
## messagesControl    0.296638   0.001058   280.4   <2e-16
## messagesCivic Duty 0.314538   0.002367   132.9   <2e-16
## messagesHawthorne  0.322375   0.002367   136.2   <2e-16
## messagesNeighbors  0.377948   0.002367   159.7   <2e-16
## 
## Residual standard error: 0.4627 on 305862 degrees of freedom
## Multiple R-squared:  0.3145, Adjusted R-squared:  0.3145 
## F-statistic: 3.508e+04 on 4 and 305862 DF,  p-value: < 2.2e-16

係数の推定値として、標本平均と同じ値が得られた。

平均因果効果は、コントロールグループとの係数の差である。

4.3.3 不均一トリートメント効果

“Neighbors”という処置について、2004年の選挙で投票した人の平均処置効果を求める。

ate_voter <- social %>%
  filter(primary2004 == 1) %>%
  filter(messages %in% c("Neighbors", "Control")) %>%
  group_by(messages) %>%
  summarize(mean_turnout = mean(primary2006)) %>%
  with(mean_turnout[2] - mean_turnout[1])
ate_voter

## [1] 0.09652525

“Neighbors”という処置について、2004年の選挙で投票しなかった人の平均処置効果を求める。

ate_nonvoter <- social %>%
  filter(primary2004 == 0) %>%
  filter(messages %in% c("Neighbors", "Control")) %>%
  group_by(messages) %>%
  summarize(mean_turnout = mean(primary2006)) %>%
  with(mean_turnout[2] - mean_turnout[1])
ate_nonvoter

## [1] 0.06929617

差を求める。

ate_voter - ate_nonvoter

## [1] 0.02722908

交差項を含む回帰モデルを推定する。

fit_int <- social %>%
  filter(messages %in% c("Neighbors", "Control")) %>%
  lm(primary2006 ~ primary2004 * messages, data = .)
summary(fit_int)

## 
## Call:
## lm(formula = primary2006 ~ primary2004 * messages, data = .)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.4823 -0.3064 -0.2371  0.6142  0.7629 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)
## (Intercept)                   0.237110   0.001345 176.322  < 2e-16
## primary2004                   0.148695   0.002125  69.963  < 2e-16
## messagesNeighbors             0.069296   0.003310  20.934  < 2e-16
## primary2004:messagesNeighbors 0.027229   0.005198   5.239 1.62e-07
## 
## Residual standard error: 0.4554 on 229440 degrees of freedom
## Multiple R-squared:  0.03078,    Adjusted R-squared:  0.03076 
## F-statistic:  2428 on 3 and 229440 DF,  p-value: < 2.2e-16

年齢変数を作成する。

social <- social %>% 
  mutate(age = 2006 - yearofbirth)
summary(social$age)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   20.00   41.00   50.00   49.79   59.00  106.00

Neighborメッセージの効果を確認するため、投票率を年齢とNeighborメッセージを受け取ったかどうかに回帰する。

fit_age <- social %>%
  filter(messages %in% c("Control", "Neighbors")) %>%
  lm(primary2006 ~ age * messages, data = .)
summary(fit_age)

## 
## Call:
## lm(formula = primary2006 ~ age * messages, data = .)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.6146 -0.3214 -0.2654  0.6227  0.8226 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)
## (Intercept)           0.0974733  0.0037603  25.922  < 2e-16
## age                   0.0039982  0.0000725  55.145  < 2e-16
## messagesNeighbors     0.0498294  0.0091519   5.445 5.19e-08
## age:messagesNeighbors 0.0006283  0.0001762   3.565 0.000364
## 
## Residual standard error: 0.4577 on 229440 degrees of freedom
## Multiple R-squared:  0.02081,    Adjusted R-squared:  0.02079 
## F-statistic:  1625 on 3 and 229440 DF,  p-value: < 2.2e-16

いくつかの年齢で、予測値を計算し、予測値のグループ毎の差として平均処置効果を求める。

age_neighbor <- data_frame(age = seq(from = 25, to = 85, by = 20),
                           messages = "Neighbors")
age_control <- data_frame(age = seq(from = 25, to = 85, by = 20),
                          messages = "Control")
## age = 25, 45, 65, 85 の平均処置効果
ate_age <- predict(fit_age, newdata = age_neighbor) - 
  predict(fit_age, newdata = age_control)
ate_age

##          1          2          3          4 
## 0.06553713 0.07810329 0.09066944 0.10323560

年齢の二乗をモデルに投入する。

fit_age2 <- social %>%
  filter(messages %in% c("Control", "Neighbors")) %>%
  lm(primary2006 ~ age + I(age^2) + messages +
       age:messages + I(age^2):messages, data = .)
summary(fit_age2)

## 
## Call:
## lm(formula = primary2006 ~ age + I(age^2) + messages + age:messages + 
##     I(age^2):messages, data = .)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.4519 -0.3344 -0.2758  0.6334  0.8749 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)
## (Intercept)                -7.385e-02  9.164e-03  -8.058 7.78e-16
## age                         1.143e-02  3.696e-04  30.914  < 2e-16
## I(age^2)                   -7.389e-05  3.605e-06 -20.494  < 2e-16
## messagesNeighbors          -4.330e-02  2.232e-02  -1.940   0.0523
## age:messagesNeighbors       4.646e-03  8.989e-04   5.169 2.36e-07
## I(age^2):messagesNeighbors -3.961e-05  8.744e-06  -4.529 5.92e-06
## 
## Residual standard error: 0.4571 on 229438 degrees of freedom
## Multiple R-squared:  0.02346,    Adjusted R-squared:  0.02344 
## F-statistic:  1102 on 5 and 229438 DF,  p-value: < 2.2e-16

年齢が25歳から85歳までについて、処置 (Neighbor) グループの予測値とコントロールグループの予測値を計算する。

## Neighborグループ
yT_hat <- data_frame(age = 25:85, messages = "Neighbors") %>% 
  predict(fit_age2, newdata = .)
## Controlグループ
yC_hat <- data_frame(age = 25:85, messages = "Control") %>% 
  predict(fit_age2, newdata = .)

結果を図示する。まずは、各グループ（各条件）について投票率の予測値を図示する。

df <- data_frame(treatment = yT_hat,
                 control = yC_hat,
                 age = 25:85) %>%
  mutate(ate = treatment - control)
pred_plt <- ggplot(df, aes(x = age)) +
  geom_line(aes(y = treatment), color = "red") +
  geom_line(aes(y = control), color = "blue", linetype = "dashed") +
  geom_text(aes(x = 50, y = 0.45, label = "Neighbors 条件"),
            color = "red", family = "HiraginoSans-W3") +
  geom_text(aes(x = 70, y = 0.325, label = "Control 条件"),
            color = "blue", family = "HiraginoSans-W3") +
  labs(x = "年齢", y = "投票率の予測値")
print(pred_plt)

次に、平均処置効果の推定値を年齢の関数として図示する。

ate_age <- ggplot(df, aes(x = age, y = ate)) +
  geom_line() +
  labs(x = "年齢", y = "平均処置効果の推定値") +
  ylim(0, 0.1)
print(ate_age)

4.3.4 回帰分断デザイン（回帰不連続デザイン）

イギリスの1950年kら1970年の間の総選挙に立候補し、勝つ見込みのあった候補者に関するデータファイル MPs.csv をダウンロードし、dataディレクトリに保存する。

download.file(url = "http://qss.princeton.press/student-files/PREDICTION/MPs.csv",
              destfile = "data/MPs.csv")

データを読み込み、中身を確認する。

MPs <- read_csv("data/MPs.csv")
glimpse(MPs)

## Observations: 427
## Variables: 10
## $ surname    <chr> "Llewellyn", "Morris", "Walker", "Walker", "Waring"...
## $ firstname  <chr> "David", "Claud", "George", "Harold", "John", "Rona...
## $ party      <chr> "tory", "labour", "tory", "labour", "tory", "labour...
## $ ln.gross   <dbl> 12.13591, 12.44809, 12.42845, 11.91845, 13.52022, 1...
## $ ln.net     <dbl> 12.135906, 12.448091, 10.349009, 12.395034, 13.5202...
## $ yob        <int> 1916, 1920, 1914, 1927, 1923, 1921, 1907, 1912, 190...
## $ yod        <int> 1992, 2000, 1999, 2003, 1989, 2002, 1987, 1984, 199...
## $ margin.pre <dbl> NA, NA, -0.057168204, -0.072508894, -0.269689620, 0...
## $ region     <chr> "Wales", "South West England", "North East England"...
## $ margin     <dbl> 0.05690404, -0.04973833, -0.04158868, 0.02329524, -...

データを二大政党別に分ける。

MPs_labour <- filter(MPs, party == "labour")
MPs_tory <- filter(MPs, party == "tory")

労働党について、選挙での勝利マージンが正の場合と負の場合に分けて回帰モデルを推定する。

labour_fit1 <- MPs_labour %>%
  filter(margin < 0) %>%
  lm(ln.net ~ margin, data = .)
labour_fit2 <- MPs_labour %>% 
  filter(margin > 0) %>% 
  lm(ln.net ~ margin, data = .)

保守党について、選挙での勝利マージンが正の場合と負の場合に分けて回帰モデルを推定する。

tory_fit1 <- MPs_tory %>%
  filter(margin < 0) %>%
  lm(ln.net ~ margin, data = .)
tory_fit2 <- MPs_tory %>% 
  filter(margin > 0) %>% 
  lm(ln.net ~ margin, data = .)

労働党の予測値を計算する。

## 負のマージンの予測範囲
y1l_range <- c(min(MPs_labour$margin), 0)
## 正のマージンの予測範囲
y2l_range <- c(0, max(MPs_labour$margin))
## 予測
y1_labour <- predict(labour_fit1,
                     newdata = data_frame(margin = y1l_range))
y2_labour <- predict(labour_fit2,
                     newdata = data_frame(margin = y2l_range))

同様に、保守党の予測値を計算する。

## 負のマージンの予測範囲
y1t_range <- c(min(MPs_tory$margin), 0)
## 正のマージンの予測範囲
y2t_range <- c(0, max(MPs_tory$margin))
## 予測
y1_tory <- predict(tory_fit1,
                   newdata = data_frame(margin = y1t_range))
y2_tory <- predict(tory_fit2,
                   newdata = data_frame(margin = y2t_range))

純資産の対数と選挙マージンの散布図に各政党の予測値を上書きして図示する。

plt_labour <- ggplot(MPs_labour) +
  geom_point(aes(x = margin, y = ln.net), pch = 16,
             color = ifelse(MPs_labour$margin > 0, "red", "blue")) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  geom_segment(x = min(MPs_labour$margin), y = y1_labour[1],
               xend = 0, yend = y1_labour[2],
               color = "blue") +
  geom_segment(x = 0, y = y2_labour[1],
               xend = max(MPs_labour$margin), yend = y1_labour[2], 
               color = "red") +
  xlim(-0.5, 0.5) + ylim(6, 18) +
  labs(x = "勝利マージン", y = "死亡時純資産の対数",
       title = "労働党")
plt_tory <- ggplot(MPs_tory) +
  geom_point(aes(x = margin, y = ln.net), pch = 16,
             color = ifelse(MPs_tory$margin > 0, "red", "blue")) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  geom_segment(x = min(MPs_tory$margin), y = y1_tory[1],
               xend = 0, yend = y1_tory[2],
               color = "blue") +
  geom_segment(x = 0, y = y2_tory[1],
               xend = max(MPs_tory$margin), yend = y1_tory[2], 
               color = "red") +
  xlim(-0.5, 0.5) + ylim(6, 18) +
  labs(x = "勝利マージン", y = "死亡時純資産の対数",
       title = "保守党")
print(plt_labour + plt_tory)

勝利マージンが0の地点における、労働党当選者と落選者の純資産の平均値を計算し、その差を求める。

(labour_won <- exp(y2_labour[1]))

##        1 
## 152603.1

(labour_lost <- exp(y1_labour[2]))

##        2 
## 220503.5

(ate_labour <- labour_won - labour_lost)

##         1 
## -67900.34

負の平均因果効果がある。

勝利マージンが0の地点における、保守党当選者と落選者の純資産の平均値を計算し、その差を求める。

(tory_won <- exp(y2_tory[1]))

##        1 
## 533813.5

(tory_lost <- exp(y1_tory[2]))

##        2 
## 278762.5

(ate_tory <- tory_won - tory_lost)

##        1 
## 255050.9

正の平均因果効果がある。

保守党候補について、プラシーボテストを実施する。

tory_fit3 <- MPs_tory %>%
  filter(margin < 0) %>%
  lm(margin.pre ~ margin, data = .)
tory_fit4 <- MPs_tory %>% 
  filter(margin > 0) %>%
  lm(margin.pre ~ margin, data = .)
coef(tory_fit4)[1] - coef(tory_fit3)[1] 

## (Intercept) 
## -0.01725578