政治学方法論 I

線形回帰 (3)

準備

いつも通り今回使うパッケージを読み込む。

library('foreign')
library('MASS')
library('dplyr')
library('ggplot2')
library('rethinking')

単回帰と重回帰

説明のために『Stataによる計量政治学』（浅野正彦, 矢内勇生. 2013）で使用されているデータ（hr96-09.dta） を使う。 まず、このデータを読み込む（データへのパスは各自の状況に応じて変えること。ここでは、RStudioのプロジェクトを利用していて、プロジェクトが存在するフォルダ内にdata という名前のフォルダがあり、その中にデータセットを保存することを仮定している）。

download.file('http://www2.kobe-u.ac.jp/~yyanai/jp/quant-methods-stata/data/hr96-09.dta',
              destfile = 'data/hr96-09.dta')
HR <- read.dta("data/hr96-09.dta")
glimpse(HR)

## Observations: 5,614
## Variables: 16
## $ year      (int) 1996, 1996, 1996, 1996, 1996, 1996, 1996, 1996, 1996...
## $ ku        (chr) "aichi", "aichi", "aichi", "aichi", "aichi", "aichi"...
## $ kun       (int) 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3...
## $ party     (fctr) NFP, LDP, DPJ, JCP, bunka, NP, msz, NFP, LDP, DPJ, ...
## $ name      (chr) "KAWAMURA, TAKASHI", "IMAEDA, NORIO", "SATO, TAISUKE...
## $ age       (int) 47, 72, 53, 43, 51, 51, 45, 51, 71, 30, 31, 44, 61, ...
## $ status    (fctr) incumbent, moto, incumbent, challenger, challenger,...
## $ nocand    (int) 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 7, 7, 7...
## $ wl        (fctr) win, lose, lose, lose, lose, lose, lose, win, lose,...
## $ rank      (int) 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3...
## $ previous  (int) 2, 3, 2, 0, 0, 0, 0, 2, 1, 1, 0, 0, 0, 0, 0, 1, 3, 1...
## $ vote      (int) 66876, 42969, 33503, 22209, 616, 566, 312, 56101, 44...
## $ voteshare (dbl) 40.0, 25.7, 20.1, 13.3, 0.4, 0.3, 0.2, 32.9, 26.4, 2...
## $ eligible  (int) 346774, 346774, 346774, 346774, 346774, 346774, 3467...
## $ turnout   (dbl) 49.22, 49.22, 49.22, 49.22, 49.22, 49.22, 49.22, 51....
## $ exp       (int) 9828097, 9311555, 9231284, 2177203, NA, NA, NA, 1294...

衆議院議員経験があることを表すダミー変数を作る。

HR <- HR %>%
    mutate(experience = as.numeric(status == "incumbent" | status == "moto"))

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

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

単回帰

1つの結果変数に対し、予測変数を1つしか使わない回帰分析を単回帰 (simple regression) と呼ぶ。以下に例を示す。

単回帰 (1): 説明変数が二値しかとらないとき（モデル1）

得票率 \(v\)（結果変数）を議員経験 \(x\)（予測変数）で説明するモデルを考えよう。 議員経験は、現職または元職の候補者なら1、そうでなければ0をとる二値 (binary) 変数（ダミー変数）である。

まず、統計モデルを作る。 \[v_i \sim \mbox{Normal}(\mu_i, \sigma)\] \[\mu_i = \alpha + \beta x_i\] \[\alpha \sim \mbox{Normal}(0, 100)\] \[\beta \sim \mbox{Normal}(0, 100)\] \[\sigma \sim \mbox{Uniform}(0, 50)\]

このモデルで推定するパラメタは、\(\alpha\), \(\beta\), \(\sigma\) の3つである。 これを rethinking::map() で推定する。

model1 <- alist(
    voteshare ~ dnorm(mu, sigma),
    mu <- alpha + beta*experience,
    alpha ~ dnorm(0, 100),
    beta ~ dnorm(0, 100),
    sigma ~ dunif(0, 50)
)
fit1 <- map(model1, data = HR09)

これで、fit1 に推定結果（係数の推定値など）が保存される。

基本的な結果は、rethinking::precis() で見ることができる。

precis(fit1)

##        Mean StdDev  5.5% 94.5%
## alpha 13.98   0.62 12.99 14.98
## beta  30.86   0.98 29.29 32.43
## sigma 16.24   0.34 15.70 16.79

この出力の、Mean の列に各パラメタのMAPの平均値が示されている。 これにより、\(\hat{\alpha} =\) 13.98, \(\hat{\beta}=\) 30.86, \(\hat{\sigma}=\) 16.24 が得られた。 したがって、 \[\hat{v}_i = 13.98 + 30.86 \cdot x_i\] と、なる。

この結果を図示しよう。推定の不確実性も同時に図示するために、まずパラメタをサンプルする。

post <- mvrnorm(10000, mu = coef(fit1), Sigma = vcov(fit1)) %>%
    as.data.frame()
head(post)

##      alpha     beta    sigma
## 1 15.16026 28.98038 16.52129
## 2 14.29087 31.79521 16.13069
## 3 13.94269 31.51089 16.78033
## 4 13.86799 31.01185 16.44082
## 5 13.96910 31.44706 16.22395
## 6 14.00164 31.51690 15.67434

パラメタのサンプルから \(\mu\) を計算する。\(x\)は0または1である。

mu <- sapply(c(0, 1), function(x) post$alpha + post$beta*x)

\(x\)の各値に対応する\(\mu\)の平均値と89%HPDIを求める。

mu_mean <- apply(mu, 2, mean)
mu_HPDI <- apply(mu, 2, HPDI)

これを図示する。

df <- data.frame(experience = c(0, 1), mu_mean,
                 lb = mu_HPDI[1, ], ub = mu_HPDI[2, ])
plt_m1 <- ggplot(HR09, aes(x = experience)) +
    geom_ribbon(data = df, aes(ymin = lb, ymax = ub), fill = 'gray') +
    geom_jitter(aes(y = voteshare), color = 'royalblue',
                position = position_jitter(width = 0.05)) +
    geom_line(data = df, aes(y = mu_mean)) +
    scale_x_continuous(breaks = c(0, 1))
print(plt_m1)

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

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

HR09 %>%
    filter(experience == 0) %>%
    with(mean(voteshare)) %>%
    round(2)

## [1] 13.98

HR09 %>%
    filter(experience == 1) %>%
    with(mean(voteshare)) %>%
    round(2)

## [1] 44.84

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

単回帰 (2): 説明変数が連続値をとるとき（モデル2, 3, 4）

同様に、得票率 \(v\) を選挙費用\(m\) (exp) で説明する統計モデルを考えよう。 モデルは次のようになる。 \[v_i \sim \mbox{Normal}(\mu_i, \sigma)\] \[\mu_i = \alpha + \gamma m_i\] \[\alpha \sim \mbox{Normal}(0, 100)\] \[\gamma \sim \mbox{Normal}(0, 100)\] \[\sigma \sim \mbox{Uniform}(0, 50)\]

このモデルで推定するパラメタは、\(\alpha\), \(\gamma\), \(\sigma\) の3つである。 これを rethinking::map() で推定する。

model2 <- alist(
    voteshare ~ dnorm(mu, sigma),
    mu <- alpha + gamma*exp,
    alpha ~ dnorm(0, 100),
    gamma ~ dnorm(0, 100),
    sigma ~ dunif(0, 50)
)
fit2 <- map(model2, data = HR09)

## Error in map(model2, data = HR09): initial value in 'vmmin' is not finite
## The start values for the parameters were invalid. This could be caused by missing values (NA) in the data or by start values outside the parameter constraints. If there are no NA values in the data, try using explicit start values.

エラーが出た。このエラーの原因は、HR09データの変数expに欠測値 (missing value) があることである。そこで、欠測値の代わりに平均値を使うことにする。

mis <- is.na(HR09$exp)
HR09_imp <- HR09
HR09_imp$exp <- ifelse(mis, mean(HR09$exp, na.rm = TRUE), HR09$exp)

欠測値が補完された新しいデータセットで、もう一度試す。

fit2 <- map(model2, data = HR09_imp)
precis(fit2)

##        Mean StdDev  5.5% 94.5%
## alpha  7.54   0.76  6.33  8.76
## gamma  0.00   0.00  0.00  0.00
## sigma 16.13   0.34 15.59 16.67

\(\gamma\)の推定値が0になっている。効果がないということだろうか？ 表示する桁数を増やしてみよう。

precis(fit2, digits = 8)

##              Mean    StdDev        5.5%       94.5%
## alpha  7.54467341 0.7587308  6.33207498  8.75727184
## gamma  0.00000307 0.0000001  0.00000292  0.00000323
## sigma 16.13001560 0.3379543 15.58989939 16.67013180

この結果をみると、\(\gamma\) は非常に小さいものの、正の効果を持っているようである。効果がこのように小さく表示されているのは、選挙費用 exp の測定単位が円であり、効果として表示されているのは、1円が得票率を何ポイント上昇させるかである。当然ながら、1円の効果はさほど大きくないだろう。こういうときは、予測変数の測定単位を変えた方が良い。試しに、測定単位を百万円に変えた exp_m という変数を使って推定し直してみよう。

HR09_imp <- HR09_imp %>%
    mutate(exp_m = exp / 10^6)
model3 <- alist(
    voteshare ~ dnorm(mu, sigma),
    mu <- alpha + gamma*exp_m,
    alpha ~ dnorm(0, 100),
    gamma ~ dnorm(0, 100),
    sigma ~ dunif(0, 50)
)
fit3 <- map(model3, data = HR09_imp)
precis(fit3)

##        Mean StdDev  5.5% 94.5%
## alpha  7.55   0.76  6.33  8.76
## gamma  3.07   0.10  2.92  3.23
## sigma 16.13   0.34 15.59 16.67

exp の代わりにexp を1,000,000 で割った exp_m を使ったので、見た目の効果量が1,000,000倍になった。もちろん、実質的な結果にはなんの影響もない。先ほどは1円の効果が示されていたのに対し、今度は100万円の効果が示されているということである。つまり、100万円支出が増えると、得票率が約3ポイント上がることが期待される。

また、\(\alpha\)の推定値である7.55 は何を意味するだろうか。これは、選挙費用が0のときの得票率の期待値である。しかし、得票率が0の候補者はいない。したがって、exp_m は分析の前に中心化したほうが良い。 exp_m を平均値で中心化して推定をやり直そう。

HR09_imp <- HR09_imp %>%
    mutate(exp_m_c = exp_m - mean(exp_m))
model4 <- alist(
    voteshare ~ dnorm(mu, sigma),
    mu <- alpha + gamma*exp_m_c,
    alpha ~ dnorm(0, 100),
    gamma ~ dnorm(0, 100),
    sigma ~ dunif(0, 50)
)
fit4 <- map(model4, data = HR09_imp)
precis(fit4)

##        Mean StdDev  5.5% 94.5%
## alpha 26.34   0.48 25.57 27.10
## gamma  3.07   0.10  2.92  3.23
## sigma 16.13   0.34 15.59 16.67

\(\gamma\)と\(\sigma\)の推定値はモデル3と同じだが、\(\alpha\) の推定値が変わった。この値は、選挙費用が平均値の候補者の得票率の期待値を示している。つまり、選挙費用について平均的な候補者の得票率は、26.34%になることが予測される。

この結果を図示しよう（図にするときは中心化していないほうがわかりやすいので、モデル3を使う）。推定の不確実性も同時に図示するために、まずパラメタをサンプルする。

post <- mvrnorm(10000, mu = coef(fit3), Sigma = vcov(fit3)) %>%
    as.data.frame()
head(post)

##      alpha    gamma    sigma
## 1 7.571149 3.052724 16.45421
## 2 7.335392 3.055441 16.47787
## 3 6.210171 3.111821 16.14939
## 4 6.973732 3.098757 15.74471
## 5 6.597751 3.129001 16.36816
## 6 7.404841 3.022677 16.11619

パラメタのサンプルから \(\mu\) を計算する。\(m\)の範囲はデータセットと同じにする。

m <- seq(min(HR09_imp$exp_m_c), max(HR09_imp$exp_m), length.out = 100)
mu <- sapply(m, function(x) post$alpha + post$gamma*x)

各m に対応する\(mu\)の平均値と89%HPDIを求める。

mu_mean <- apply(mu, 2, mean)
mu_HPDI <- apply(mu, 2, HPDI)

これを図示する。

df <- data.frame(exp_m = m, mu_mean,
                 lb = mu_HPDI[1, ], ub = mu_HPDI[2, ])
plt_m3 <- ggplot(HR09_imp, aes(x = exp_m)) +
    geom_ribbon(data = df, aes(ymin = lb, ymax = ub), fill = 'gray') +
    geom_jitter(aes(y = voteshare), color = 'royalblue',
                position = position_jitter(width = 0.05)) +
    geom_line(data = df, aes(y = mu_mean)) +
    xlab('expenditure (million yen)')
print(plt_m3)

重回帰

上の例を見ると、議員経験と選挙費用の両者が得票率に影響しているようである。どちらも結果変数の値の変化を説明していると思われるので、予測変数として2つの変数を使ってみよう。以下で行うような予測変数が2つ以上の回帰分析は、重回帰 (multiple regression) と呼ばれる。

重回帰の係数の推定

まず、統計モデルを作ろう。以下のような線形回帰モデルを考えることにする。 \[v_i \sim \mbox{Normal}(\mu_i, \sigma)\] \[\mu_i = \alpha + \beta x_i + \gamma m_i\] \[\alpha \sim \mbox{Normal}(0, 100)\] \[\beta \sim \mbox{Normal}(0, 100)\] \[\gamma \sim \mbox{Normal}(0, 100)\] \[\sigma \sim \mbox{Uniform}(0, 50)\]

推定するパラメタは、\(\alpha\), \(\beta\), \(\gamma\), \(\sigma\) の4つである。 これをmap() で推定しよう。

model5 <- alist(
    voteshare ~ dnorm(mu, sigma),
    mu <- alpha + beta*experience + gamma*exp_m_c,
    alpha ~ dnorm(0, 100),
    beta ~ dnorm(0, 100),
    gamma ~ dnorm(0, 100),
    sigma ~ dunif(0, 50)
)
fit5 <- map(model5, data = HR09_imp)
precis(fit5)

##        Mean StdDev  5.5% 94.5%
## alpha 19.12   0.66 18.07 20.17
## beta  18.02   1.22 16.06 19.98
## gamma  1.86   0.12  1.66  2.05
## sigma 14.78   0.31 14.29 15.28

結果をキャタピラプロットで示したいときは、precis() をplot()すればよい（ただい、論文に載せるにはもう少し改善が必要である）。

plot(precis(fit5))

この結果の意味を考えてみよう。まず、切片であるalpha の推定値19.13は、experience と exp_m_cの両者が0のときの予測得票率を示す。言い換えると、選挙費用が平均値に一致する新人候補の予測得票率は19.13% である。ただし、選挙費用が平均値に一致する新人候補が実際に存在するとは限らないことに注意する必要がある。

次に、betaの値の意味を考えてみよう。この値は、他の変数を一定にしたとき、議員経験のある候補者と議員経験がない候補者の得票率の差を表す。つまり、選挙費用が同じ候補者を比べれば、経験のある候補者の得票率は新人候補の得票率より約18ポイント高いことが期待される。ただし、選挙費用がまったく同じ新人候補と現職あるいは元職候補が存在するとは限らない。

最後に、gammaの意味を考える。この値は、他の変数を一定にしたとき、選挙費用が1単位増えるごとに得票率がどの程度上がるかを示している。つまり、議員経験が同じ候補者を比べれば、選挙費用が100万円増えるごとに得票率が1.86ポイント上昇することが期待される。

このように、重回帰の係数は、他の予測変数を一定にしたときに該当する予測変数が結果変数に与える影響を示している。

論文に結果を書くときは、それぞれの効果を実質的に議論する必要がある。例えば、betaの推定値は18なので、新人候補はかなり不利な戦いを強いられているということであり、この効果は実質的に意味があるものと言えるだろう。

重回帰の結果を図示する

単回帰では、結果変数1つと予測変数1つの関係を分析するので、結果変数を縦軸、予測変数を横軸にした二次元の図で結果を示すことが容易にできた。 しかし、重回帰では予測変数が2つ以上あるので、少し工夫が必要になる。ここでは、3つの方法を紹介する。

(1) 予測変数の残差を利用した図示

重回帰を行う1つの理由は、ある予測変数では説明しきれない部分を他の予測変数で説明しようとすることである（他の理由については次回または方法論 II で解説する予定）。 結果変数 \(v\) を説明するために\(x\)だけでなく\(m\) も加えることで、

1. すでに\(x\)の値は知っている状況で、\(m\)の値もわかったらどれだけ\(v\)の予測がよくなるか、
2. すでに\(m\)の値は知っている状況で、\(x\)の値もわかったらどれだけ\(v\)の予測がよくなるか

を推定することができる。つまり、各変数の影響は、他の変数の影響を除去した後に残るものである。よって、他の変数の影響を除去してみよう。

まず、議員経験に注目し、選挙費用の影響を除去する。そのために、まず議員経験を選挙費用に回帰し、議員経験に含まれる選挙費用の影響を取り除く。 以下の統計モデルを使う。 \[x_i \sim \mbox{Normal}(\mu_i, \sigma)\] \[\mu_i = a_1 + b_1 m\] \[a_1 \sim \mbox{Normal}(0, 100)\] \[b_1 \sim \mbox{Normal}(0, 100)\] \[\sigma \sim \mbox{Uniform}(0, 50)\]

\(a_1\), \(b_1\), \(\sigma\) を推定する。

res_model1_1 <- alist(
    experience ~ dnorm(mu, sigma),
    mu <- a1 + b1*exp_m_c,
    a1 ~ dnorm(0, 100),
    b1 ~ dnorm(0, 100),
    sigma ~ dunif(0, 50)
)
res_fit1_1 <- map(res_model1_1, data = HR09_imp)

この推定結果を用い、各候補者の得票率の予測値 (\(\mu\)) を計算する。

mu <- coef(res_fit1_1)[1] + coef(res_fit1_1)[2]*HR09_imp$exp_m_c

この予測値を使い、残差 (residual) を計算する。残差は、観測値と予測値の差である。

HR09_imp$experience_res <- HR09_imp$experience - mu

次に、結果変数である \(v\) をこの残差に回帰する。 以下の統計モデルを使おう。 \[v_i \sim \mbox{Normal}(x_i^{\mathrm{res}}, \sigma)\] \[\mu_i = \rho + \beta x_i^{\mathrm{res}}\] \[\rho \sim \mbox{Normal}(0, 100)\] \[\beta \sim \mbox{Normal}(0, 100)\] \[\sigma \sim \mbox{Uniform}(0, 50)\]

パラメタを推定する。

res_model1_2 <- alist(
    voteshare ~ dnorm(mu, sigma),
    mu <- rho + beta*experience_res,
    rho ~ dnorm(0, 100),
    beta ~ dnorm(0, 100),
    sigma ~ dunif(0, 50)
)
res_fit1_2 <- map(res_model1_2, data = HR09_imp)
precis(res_fit1_2)

##        Mean StdDev  5.5% 94.5%
## rho   26.34   0.63 25.33 27.34
## beta  18.01   1.76 15.20 20.82
## sigma 21.23   0.44 20.52 21.95

このbetaの推定値は、重回帰で求めた推定値と（計算誤差を除いて）一致する！これは偶然ではない。上で述べた通り、重回帰は他の変数の影響を除去しているので、このように残差を使った回帰でも同じ結果にいたる。

この2段階目の回帰分析は単回帰なので、単回帰の場合と同じように結果を図示することができる。

post <- mvrnorm(10000, mu = coef(res_fit1_2), Sigma = vcov(res_fit1_2)) %>%
    as.data.frame()
x <- seq(min(HR09_imp$experience_res), max(HR09_imp$experience_res), 
         length.out = 100)
mu <- sapply(x, function(x) post$rho + post$beta*x)
mu_mean <- apply(mu, 2, mean)
mu_HPDI <- apply(mu, 2, HPDI)  ## 89% HPDI
df <- data.frame(experience_res = x, mu_mean,
                 lb = mu_HPDI[1,], ub = mu_HPDI[2,])
plt_res1 <- ggplot(HR09_imp, aes(x = experience_res)) +
    geom_ribbon(data = df, aes(ymin = lb, ymax = ub), fill = 'gray') +
    geom_point(aes(y = voteshare), color = 'royalblue') +
    geom_line(data = df, aes(y = mu_mean)) +
    xlab('Experinece residuals') +
    geom_vline(xintercept = 0, linetype = 'dashed')
print(plt_res1)

上の図で、垂直な点線は残差0の位置を示している。この点線より左側の点は、選挙費用で予測されるよりも議員経験が少ない観測値である。また、右側にあるのは、選挙費用から予測されるよりも議員経験がある候補者の点である。この図から、議員経験が得票率に正の影響を及ぼしていることがわかる。

次に、選挙費用の残差を横軸にして、これと同様の図を作ってみよう（今週の課題を参照）。

(2) 予測変数の周辺効果の図示

次に、他の変数を特定の値に固定して、1つの予測変数が結果変数に与える影響、つまり、特定の予測変数の周辺効果 (marginal effect) を図示してみよう。

ここでは、議員経験の周辺効果を考える。 議員経験が取り得る値は0か1だが、図を作るために0と1の間の値についても予測値を計算することにする。

df <- data.frame(experience = seq(0, 1, length.out = 100))

次に、選挙費用をどの値に固定するか決める。とりあえず、平均値に固定してみよう。

df$money_mean <- with(HR09_imp, mean(exp_m_c))  ## this must be 0

これらの値を使って、予測値 \(\mu\) を計算する。

post <- mvrnorm(10000, mu = coef(fit5), Sigma = vcov(fit5)) %>%
    as.data.frame()
mu <- sapply(df$experience,
             function(x) post$alpha + post$beta*x + post$gamma*df$money_mean)
# Because money_mean = 0, you can run the following instead.
# mu <- sapply(df$x, function(x) post$alpha + post$beta*x)
df$mu_mean <- apply(mu, 2, mean)
mu_HPDI <- apply(mu, 2, HPDI)  # 89% HPDI
df$lb <- mu_HPDI[1,]
df$ub <- mu_HPDI[2,]
plt_x <- ggplot(HR09_imp, aes(x = experience)) +
    geom_ribbon(data = df, aes(ymin = lb, ymax = ub), fill = 'gray') +
    geom_jitter(aes(y = voteshare), color = 'royalblue',
                position = position_jitter(width = 0.05)) +
    geom_line(data = df, aes(y = mu_mean)) +
    scale_x_continuous(breaks = c(0, 1))
print(plt_x)

他の変数は、平均値意外に設定しても良い。 例えば、最大値と最小値も使って、次ようなことができる。

df$money_min <- with(HR09_imp, min(exp_m_c))
df$money_max <- with(HR09_imp, max(exp_m_c))
mu_min <- sapply(df$experience,
                 function(x) post$alpha + post$beta*x + post$gamma*df$money_min)
mu_max <- sapply(df$experience,
                 function(x) post$alpha + post$beta*x + post$gamma*df$money_max)
df$mu_min_mean <- apply(mu_min, 2, mean)
df$mu_max_mean <- apply(mu_max, 2, mean)
mu_min_HPDI <- apply(mu_min, 2, HPDI)  # 89% HPDI
mu_max_HPDI <- apply(mu_max, 2, HPDI)  # 89% HPDI
df$lb_min <- mu_min_HPDI[1,]
df$ub_min <- mu_min_HPDI[2,]
df$lb_max <- mu_max_HPDI[1,]
df$ub_max <- mu_max_HPDI[2,]
plt_x <- ggplot(HR09_imp, aes(x = experience)) +
    geom_ribbon(data = df, aes(ymin = lb, ymax = ub), fill = 'gray') +
    geom_ribbon(data = df, aes(ymin = lb_min, ymax = ub_min), fill = 'gray') +
    geom_ribbon(data = df, aes(ymin = lb_max, ymax = ub_max), fill = 'gray') +
    geom_jitter(aes(y = voteshare), position = position_jitter(width = 0.05)) +
    geom_line(data = df, aes(y = mu_mean)) +
    geom_line(data = df, aes(y = mu_min_mean), color = 'blue') +
    geom_line(data = df, aes(y = mu_max_mean), color = 'red') +
    scale_x_continuous(breaks = c(0, 1))
print(plt_x)

この図で、黒い直線は選挙費用が平均値の場合、赤は選挙費用が最大値の場合、青は選挙費用が最小値の場合の予測得票率を示している。

同じように、選挙費用の周辺効果を図にしてみよう（今週の課題を参照）。

(3) 事後予測を使った図示

最後に、この統計モデルが観察されたデータにフィットするかどうかを確かめるため、結果変数の事後予測値 (posterior predictions) と観測値の散布図を作る。

まず、予測変数の観測値を使って結果変数の事後予測値を計算する。

predictor <- list()
for(i in 1:nrow(HR09_imp)) {
    row_i <- HR09_imp[i,]
    predictor[[i]] <-  c(row_i$experience, row_i$exp_m_c)
}
mu <- sapply(predictor,
             function(x) post$alpha + post$beta*x[1] + post$gamma*x[2])
HR09_imp$mu_mean <- apply(mu, 2, mean)
mu_HPDI <- apply(mu, 2, HPDI)  # 89% HPDI
HR09_imp$mu_lb <- mu_HPDI[1,]
HR09_imp$mu_ub <- mu_HPDI[2,]

これを図示する。

pred_obs <- ggplot(HR09_imp, aes(voteshare, mu_mean)) +
    geom_point(size = .5) +
    geom_pointrange(aes(ymin = mu_lb, ymax = mu_ub)) +
    xlim(0, 100) + ylim(0, 100) +
    geom_abline(intercept = 0, slope = 1, linetype = 'dashed') +
    labs(x = 'observed vote share', y = 'predicted vote share')
print(pred_obs)

多項式を使った線形回帰

予測変数と結果変数の関係が、常に直線的であるとは限らない。 例えば、前回も使った身長のデータで、体重と身長の関係は直線的とは言えない。図で確認してみよう。

data(Howell1)
myd <- Howell1
scat <- ggplot(myd, aes(weight, height)) + 
    geom_point() +
    geom_smooth()
print(scat)

この図からわかるように、成人になるまでの身長は体重とともに伸びる。しかし、成人なると成長が止まり、体重が増えても身長はそれほど増えない（身長の成長がともっても体重だけは増やせる！）。

このように、ある予測変数（年齢）が結果変数（身長）に与える影響は、線形とは限らない。線形回帰は、予測変数と結果変数の関係が線形関数で表されることを想定しているが、多項式を使えば、曲線的な関係も表現することができる。例えば、二次の多項式を使って、身長 (\(h\)) と体重 (\(x\)) の関係を次のような統計モデルで表すことができる。 \[h_i \sim \mbox{Normal}(\mu_i, \sigma)\] \[\mu_i = \alpha + \beta_1 x_i + \beta_2 x_i^2 \] \[\alpha \sim \mbox{Normal}(168, 100)\] \[\beta_1 \sim \mbox{Normal}(0, 10)\] \[\beta_2 \sim \mbox{Normal}(0, 10)\] \[\sigma \sim \mbox{Uniform}(0, 50)\]

このモデルで、体重 \(x\) は、\(x\) としてだけでなく、\(x^2\)としても線形モデルに登場する。 \(\mu\) は \(x\) の二次関数なので、曲線的な関係が表現できる。一方、\(z = x^2\) とすれば、 \[\mu_i = \alpha + \beta_1 x_i + \beta_2 z_i\] なので、線形回帰問題と考えることが可能である（ただし、後で述べるように解釈はややこしくなる）。

このモデルを使ってパラメタを推定してみよう。

# 予測変数の準備
myd <- myd %>%
    mutate(weight_c = weight - mean(weight),  # centering
           weight_c_sq = weight_c^2)          # squared term
model5 <- alist(
    height ~ dnorm(mu, sigma),
    mu <- alpha + beta1*weight_c + beta2*weight_c_sq,
    alpha ~ dnorm(168, 100),
    beta1 ~ dnorm(0, 10),
    beta2 ~ dnorm(0, 10),
    sigma ~ dunif(0, 50)
)
fit5 <- map(model5, data = myd)

推定結果を見てみよう。

precis(fit5)

##         Mean StdDev   5.5%  94.5%
## alpha 146.66   0.37 146.06 147.26
## beta1   1.45   0.02   1.42   1.49
## beta2  -0.04   0.00  -0.04  -0.04
## sigma   5.75   0.17   5.47   6.03

ここで、それぞれの係数は何を意味しているだろうか。 予測変数を中心化したので、切片である alpha の値は、体重が平均のときの身長の平均値を表している。

では、beta1 や beta2 の値は何を意味しているだろうか。 beta1 の平均値である1.45 は、他の変数が一定のとき体重が1kg増えると身長が1.45cm増えることを意味しているはずだが、そのようなことは現実にあるだろうか？もちろん、そのようなことはありえない。ここで「他の変数」というのは「体重の二乗」のことであり、体重の二乗の値を一定に保ったまま体重自体を1kg増やすことは不可能である。したがって、多項式（他の変数の二乗）を使った線形回帰では、結果を表で提示することはいつも以上に不親切である。表から結果が意味する実質的な効果を読み取ることは非常に困難である。

よって、（多項式でない場合もそうだが、多項式の場合は特に）結果を図示することが重要になる。

まず、パラメタの値を100個シミュレートする。

post <- mvrnorm(10000, mu = coef(fit5), Sigma = vcov(fit5)) %>%
    as.data.frame()

次に、シミュレートされた値を用い、各体重に対応する\(\mu\)を計算する。

w_vec <- seq(min(myd$weight_c), max(myd$weight_c), length.out = 100)
mu_w_w2 <- sapply(w_vec,
                  function(x) post$alpha + post$beta1*x + post$beta2*(x^2))

得られた\(\mu\)の平均値（回帰曲線）と93%HPDIを求める。

mu_mean <- apply(mu_w_w2, 2, mean)
mu_HPDI <- apply(mu_w_w2, 2, HPDI, prob = .93)  # 93% HPDI

回帰曲線と推定された曲線の不確実性を図示する。中心化した値よりも元の値の方がわかりやすいので、横軸を scale_x_continuous() で修正する。

df <- data.frame(weight_c = w_vec, mu_mean,
                 lb = mu_HPDI[1, ], ub = mu_HPDI[2, ])
range(myd$weight)

## [1]  4.252425 62.992589

intrvl <- ggplot(myd, aes(x = weight_c)) +
    geom_ribbon(data = df, aes(ymin = lb, ymax = ub), fill = 'gray60') +
    geom_point(aes(y = height), color = 'royalblue') +
    geom_line(data = df, aes(y = mu_mean)) +
    scale_x_continuous(name = 'weight',
                       breaks = seq(0, 60, by = 10) - mean(myd$weight),
                       labels = seq(0, 60, by = 10)) +
    labs(x = 'weight (kg)', y = 'height (cm)')
print(intrvl)

このように、体重と身長の非直線的な関係を図示できる。 この図から、だいたい50kgを超えると、体重と身長の正の関係ががなくなることがわかる。推定された二次式を変形すると、 \[\mu = 146.66 + 1.45x - 0.04x^2 \approx -0.04 (x - 18.13)^2 + 159.81\] なので、中心化された体重が約18.13kgのときに効果が正から負に変わる。 平均体重は35.61kg だから、これを足せば、53.74kgで効果の向き（符号）に変化があることがわかる。

授業の内容に戻る