第14章　ロジスティック回帰分析

14.0 準備

準備として、作業ディレクトリを指定する。作業ディレクトリの指定にはsetwd()関数を利用する。ディレクトリへのパス（引用符の中身）は、各自で設定する。

setwd("~/projects/Stata-R/")

この章で使うパッケージをあらかじめ読み込んでおく。

library("foreign")
library("dplyr")
library("ggplot2")
## Windows ユーザは次の行をコメントアウト
theme_set(theme_gray(base_size=12, base_family="HiraKakuProN-W3"))
## 必要なら、reshape2をインストール
# install.packages('reshape2', dependencies = TRUE)
library("reshape2")

14.1 ロジスティック関数

教科書 p.250 にあるとおり、ロジスティック関数とは、何かが起きる (\(y = 1\)) か起きないか (\(y = 0\)) を考えているときに、 \[\mathrm{Pr}(y_i=1) = p_i = \mathrm{logistic}(x_i) = \frac{\exp(x_i)}{1 + \exp(x_i)} = \frac{1}{1 + \exp(-x_i)}\] となる関数である。ちなみに、ロジット関数は \[x_i = \mathrm{logit}(p_i) = \log \left( \frac{p_i}{1-p_i}\right)\] であり、ロジスティック関数の逆関数である。統計学の教科書等では、\(\mathrm{logistic}(x)\) ではなく、\(\mathrm{logit}^{-1}(x)\)と書くのが普通である。

まず、ロジスティック関数（ロジット関数の逆関数）をRで定義してみよう。 function() を使えば、自分で関数を作ることができる。

invlogit <- function(x){ ## ロジスティック関数を定義する
    ## 引数：x = 数値ベクトル
    ## 返り値：y
    return( y <- 1 / (1 + exp(-x)) )
}

これでinvlogit() 関数（ロジスティック関数）ができた。

確認のため、\(x \in [-5, 5]\) の範囲で\(\mathrm{logistic}(x) = \mathrm{logit}^{-1}(x)\)の値を図示してみよう。

x <- seq(-5, 5, length = 100)
plot(x, invlogit(x), type = "l", lwd = 2, col = "darkblue",
     ylab = "invlogit(x)", main = "ロジスティック関数")
abline(h = c(0, 1), col = "gray")

14.2 ロジスティック回帰分析の手順

例14-1 (p.252) 同様、候補者の当選 (\(y=1\)) と落選 (\(y=0\)) の違いを選挙費用と当選回数で説明するモデルを考える。ただし、ここで使用するデータ（ml_wl.dta）は架空のものである。

14.2.1 帰無仮説と対立仮説を設定する

ここで考える帰無仮説は、

\(\mathrm{H}_{01}\)：当選回数は小選挙区での当選確率に影響を与えない
\(\mathrm{H}_{02}\)：選挙費用は小選挙区での当選確率に影響を与えない

という2つの仮説である。また、対立仮説は、

\(\mathrm{H}_{a1}\)：当選回数が多いほど、小選挙区での当選確率が高くなる
\(\mathrm{H}_{a2}\)：選挙費用が多いほど、小選挙区での当選確率が高くなる

という2つの仮説である。

14.2.2 説明変数と応答変数の散布図

まず、データ（ml_wl.dta）を読み込む。

ML.WL <- read.dta("data/ml_wl.dta")

このデータを使い、自民党候補の小選挙区における当落 (wlsmd) を縦軸に、それまでの当選回数 (previous) と選挙費用 (expm) を横軸にした散布図を表示してみよう。

p1 <- ggplot(ML.WL, aes(previous, wlsmd)) + geom_point()
p2 <- ggplot(ML.WL, aes(expm, wlsmd)) + geom_point()
print(p1)

print(p2)

これで図14.3, 14.4 (p.255) と同じような図ができた。

また、相関係数は、

with(ML.WL, cor(previous, wlsmd))

## [1] 0.3500396

with(ML.WL, cor(expm, wlsmd))

## [1] 0.6132906

である。

14.2.3 ロジスティック回帰式の推定値を求める

ロジスティック回帰分析を行うにはglm() 関数を利用する。 glm() 関数はロジスティク回帰以外にも様々な一般化線形モデルで利用される。ロジスティック回帰に使うには、family を以下のように指定する。

fit.wl1 <- glm(wlsmd ~ previous + expm, data = ML.WL,
               family = binomial(link = "logit"))
summary(fit.wl1)

## 
## Call:
## glm(formula = wlsmd ~ previous + expm, family = binomial(link = "logit"), 
##     data = ML.WL)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5741  -0.3781   0.2013   0.3943   1.4948  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  -6.3811     3.5147  -1.816   0.0694
## previous      0.8085     0.5851   1.382   0.1670
## expm          0.8088     0.4000   2.022   0.0431
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 20.728  on 14  degrees of freedom
## Residual deviance: 10.384  on 12  degrees of freedom
## AIC: 16.384
## 
## Number of Fisher Scoring iterations: 6

これで、表14.2 (p.257) と同じ結果が得られた。

14.2.4 ロジスティック回帰式の精度を求める

教科書同様、当選確率の予測値 \(\hat{p}\) を求めてみよう。 Rでは、通常はpredict() 関数を使うと、予測値を計算してくれる。この関数を使ってみよう。

predicted.wl1 <- predict(fit.wl1)
## 代わりに次のようにしても同じである
## predicted.wl1 <- fit.wl1$linear.predictor
predicted.wl1

##           1           2           3           4           5           6 
##  1.70706688  2.51559946  3.24296469  3.88944860  3.64623297 -0.72052766 
##           7           8           9          10          11          12 
## -0.72081292  2.51474369 -4.60171222 -2.90319464  0.89682276 -2.33759282 
##          13          14          15 
##  0.08829018 -3.79317964 -0.72024240

これで確率の予測値が求められたのだろうか。上に表示された数字を見れば分かるとおり、これは確率ではない（0より小さい数字や1より大きい数字がある）。では、これらの数字は何を表しているのだろうか。

ここで、fit.wl1 をどうやって推定したかを思い出してみよう。推定にはglm() 関数を利用したが、そのとき link = "logit" を指定した。簡単にいうと、link = "logit" を指定したということは、この回帰式はロジット関数で表現されているということである。そして、上に表示されているpredicted.wl1 という予測値は、ロジットの値である。

これを確率に直すには、既に作ったinvlogit() 関数を使えばよい。

(predicted.wl1 <- invlogit(predicted.wl1))

##           1           2           3           4           5           6 
## 0.846455458 0.925228190 0.962419484 0.979953462 0.974574119 0.327276799 
##           7           8           9          10          11          12 
## 0.327213998 0.925168966 0.009934946 0.051995866 0.710296141 0.088057027 
##          13          14          15 
## 0.522058217 0.022027721 0.327339606

これらの値は、図14.5 (p.258) の値と一致する。

predict() 関数を使って同じ値を求めるには、type = "response" を指定すればよい。

predict(fit.wl1, type = "response")

##           1           2           3           4           5           6 
## 0.846455458 0.925228190 0.962419484 0.979953462 0.974574119 0.327276799 
##           7           8           9          10          11          12 
## 0.327213998 0.925168966 0.009934946 0.051995866 0.710296141 0.088057027 
##          13          14          15 
## 0.522058217 0.022027721 0.327339606

これは、invlogit()関数を使って求めた値と一致する。

次に、予測の正解率を求めてみよう。先ほど求めた予測値を使って、次のように予測の正解率を出すことができる。

error.rate.wl1 <- mean((predicted.wl1 > 0.5 & ML.WL$wlsmd == 0) | 
                           predicted.wl1 < 0.5 & ML.WL$wlsmd ==1)
1 - error.rate.wl1

## [1] 0.7333333

これで、表14.3 (p.259) と同様の結果が得られた。

14.2.5 回帰係数を包括的に検定する

表14.2 (p.257) にある“LR chi2(2)” に相当する数字は、Rの分析結果に表示されていない。この数字は、summary(fit.wl1) で表示される結果に示される Null deviance と Residual deviance の差である。また、“LR chi2(2)” の括弧の中にある2は自由度であるが、これはNull deviance の自由度とResidual deviance の自由度の差である。これらの数字を表から読み取って計算しても良いが、次のようにすると、これらの数値を求めることができる。

chisq.wl1 <- with(fit.wl1, null.deviance - deviance)
df.wl1 <- with(fit.wl1, df.null - df.residual)

これらの数値を使って\(\chi^2\)検定の\(p\)値を求めよう。

pchisq(chisq.wl1, df.wl1, lower.tail = FALSE)

## [1] 0.005674637

これで、表14.2 と同じ\(p\)値 (0.0057) が得られた。

14.2.6 回帰係数を個別に検定する

それぞれの変数の影響が0といえるかどうかを個別に検定するには、係数の\(p\)値を利用する。 Rでは以下の結果に現れる“Pr(>|z|)” の欄の値が係数の\(p\)値である。

summary(fit.wl1)

## 
## Call:
## glm(formula = wlsmd ~ previous + expm, family = binomial(link = "logit"), 
##     data = ML.WL)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5741  -0.3781   0.2013   0.3943   1.4948  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  -6.3811     3.5147  -1.816   0.0694
## previous      0.8085     0.5851   1.382   0.1670
## expm          0.8088     0.4000   2.022   0.0431
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 20.728  on 14  degrees of freedom
## Residual deviance: 10.384  on 12  degrees of freedom
## AIC: 16.384
## 
## Number of Fisher Scoring iterations: 6

したがって、previous（当選回数）と expm（選挙費用）の\(p\)値はそれぞれ0.167, 0.043 であり、教科書p.262 と同じ結果が示された。

14.2.7 説明変数を複数の値に固定したグラフを作成する

まず、表14.4 (p.264) と同じ表を作ってみよう。 previous と expm のそれぞれについて、調べる対象となる値のベクトルを作る。必要な値は、

最小値
最大値
0
1
平均値 - 1/2
平均値 + 1/2
平均値 - sd/2
平均値 + sd/2

の8つの値である。

prev.mean <- mean(ML.WL$previous)  # previousの平均値
prev.sd <- sd(ML.WL$previous)      # previousのsd
# previousの8つの値
prev.vec <- c(min(ML.WL$previous), max(ML.WL$previous), 0, 1,
              prev.mean - 1/2, prev.mean + 1/2,
              prev.mean - prev.sd/2, prev.mean + prev.sd/2)
expm.mean <- mean(ML.WL$expm)  # expmnの平均値
expm.sd <- sd(ML.WL$expm)      # expmのsd
## expmの8つの値
expm.vec <- c(min(ML.WL$expm), max(ML.WL$expm), 0, 1,
              expm.mean - 1/2, expm.mean + 1/2,
              expm.mean - expm.sd/2, expm.mean + expm.sd/2)
k <- 8  # 調べる値の数

次に、previous, expm それぞれについて、他方を平均値に固定した上で、8つの値を取ったときの当選確率の予測値を計算する。

## previousの値を変えて当選確率を予測する：expmは平均値に固定
wl.prev <- predict(fit.wl1, type = "response",
                   newdata = list( previous = prev.vec, expm = rep(expm.mean, k)))
## expmの値を変えて当選確率を予測する：previousは平均値に固定
wl.expm <- predict(fit.wl1, type = "response",
                   newdata = list(previous = rep(prev.mean, k), expm = expm.vec))
## 結果をまとめる
res <- rbind(wl.prev, wl.expm)

最後に、変化量を調べるために差をとる。たとえば、previous が最小値から最大値に変化したときの当選確率の変化を見るために、previousが最大値のときの予測値から、最小値のときの予測値を引く。

tbl.14.4 <- cbind(res[,2] - res[,1],
                  res[,4] - res[,3],
                  res[,6] - res[,5],
                  res[,8] - res[,7])
colnames(tbl.14.4) <- c("min -> max", "0 -> 1", "-+1/2", "-+sd/2")
tbl.14.4

##         min -> max     0 -> 1     -+1/2    -+sd/2
## wl.prev  0.8692498 0.10002099 0.1978768 0.4124971
## wl.expm  0.8995518 0.02132257 0.1979448 0.5143434

これで、表14.4 (p.264) のうち、私たちが知りたい情報を得ることができた。

当選0回の候補者に関する条件付き効果プロット

当選回数 (previous) が0のとき、選挙費用 (expm) の値の変化によって当選確率がどのように変化するかを示す図を作る。

まず、選挙費用の変化に応じた予測値を求める。求める範囲は、選挙費用 (expm) の観測された最小値と最大値の間とする。

ML.WL <- mutate(ML.WL,
                exp.vec = seq(min(expm), max(expm), length = n()),
                p.hat.wl.prev0 = predict(fit.wl1, type = "response",
                                         list(previous = rep(0, n()), expm = exp.vec)))

これを図示してみよう。

wl.exp.lab <- labs(x = "選挙費用（10万円）", y = "当選確率")
p.hat.0 <- ggplot(ML.WL, aes(exp.vec, p.hat.wl.prev0)) + geom_line(color = "red")
p.hat.0 + wl.exp.lab + ggtitle("当選回数が0の場合")

ggplot2 を使わない場合は、

plot(p.hat.wl.prev0 ~ exp.vec, data = ML.WL,
     type = "l", lwd = 2, col = "violetred",
     xlab = "選挙費用（100万円）", ylab = "当選確率",
     main = "当選回数が0の場合")

これで図14.7 (p.265) と同じような図ができた。

当選5回の候補者に関する条件付き効果プロット

当選回数 (previous) が5のとき、選挙費用 (expm) の値の変化によって当選確率がどのように変化するかを示す図を作る。

まず、選挙費用の変化に応じた予測値を求める。求める範囲は、選挙費用 (expm) の観測された最小値と最大値の間とする。

ML.WL <- mutate(ML.WL,
                p.hat.wl.prev5 = predict(fit.wl1, type = "response",
                                         list(previous = rep(5, n()), expm = exp.vec)))

これを図示してみよう。

p.hat.5 <- ggplot(ML.WL, aes(exp.vec, p.hat.wl.prev5)) + geom_line(color = "blue")
p.hat.5 + wl.exp.lab + ggtitle("当選回数が5の場合")

ggplot2を使わない場合は、

plot(p.hat.wl.prev5 ~ exp.vec, data = ML.WL,
     type = "l", lwd = 2, col = "darkblue",
     xlab = "選挙費用（100万円）", ylab = "当選確率",
     main = "当選回数が5の場合")

これで図14.8 (p.267) と同じような図ができた。

これら2つの図をひとつにまとめよう。まず、ggplot2で使いやすい形のデータフレームP.hat.wl を作る。

P.hat.wl <- melt(ML.WL, id = "exp.vec", 
                  measure = c("p.hat.wl.prev0", "p.hat.wl.prev5"))
head(P.hat.wl)

##    exp.vec       variable       value
## 1 2.000000 p.hat.wl.prev0 0.008463644
## 2 2.571429 p.hat.wl.prev0 0.013369835
## 3 3.142857 p.hat.wl.prev0 0.021059634
## 4 3.714286 p.hat.wl.prev0 0.033024260
## 5 4.285714 p.hat.wl.prev0 0.051429248
## 6 4.857143 p.hat.wl.prev0 0.079250980

tail(P.hat.wl)

##      exp.vec       variable     value
## 25  7.142857 p.hat.wl.prev5 0.9688960
## 26  7.714286 p.hat.wl.prev5 0.9801792
## 27  8.285714 p.hat.wl.prev5 0.9874225
## 28  8.857143 p.hat.wl.prev5 0.9920403
## 29  9.428571 p.hat.wl.prev5 0.9949713
## 30 10.000000 p.hat.wl.prev5 0.9968265

このデータフレームを利用し、図を作ろう。

p.hat.wl <- ggplot(P.hat.wl, aes(exp.vec, value)) + wl.exp.lab
p.hat.wl <- p.hat.wl + geom_line(aes(color = variable))
p.hat.wl + scale_color_discrete(name="当選回数", labels=c("0回", "5回")) +
    guides(color = guide_legend(reverse = TRUE))

ggplot2を使わない場合は、

plot(p.hat.wl.prev0 ~ exp.vec, data = ML.WL,
     type = "l", lwd = 2, col = "violetred", ylim = c(0, 1),
     xlab = "選挙費用（100万円）", ylab = "当選確率", main = "")
lines(p.hat.wl.prev5 ~ exp.vec, data = ML.WL, lwd = 2, col = "darkblue")
abline(h = c(0, 1), col = "gray")
legend("bottomright", title = "当選回数", legend = c("5回", "0回"),
       col = c("darkblue", "violetred"), lwd = 2, bg = "white")

これで図14.9 (p.268) と同じような図ができた。

14.3 衆院選データを使ったロジスティック回帰分析

今度は、現実の衆院選データ（hr96-09.dta）を使い、2009年の衆院選での自民党候補の小選挙区での当落を、選挙費用と当選回数で説明することを試みる。

14.3.1 帰無仮説と対立仮説

Rを使う必要がないので省略。

14.3.2 説明変数と応答変数の散布図

まず、使用するデータを用意する。

HR <- read.dta("data/hr96-09.dta")
HR <- mutate(HR,
             wlsmd = as.numeric(wl == "win"),
             expm = exp / 10^6)
HR09LDP <- filter(HR, year == 2009, party == "LDP")
## expが欠測の個体が1つあるので、それを除外しておく
HR09LDP <- filter(HR09LDP, !is.na(exp))

次に説明変数と応答変数の関係を示す散布図を作る。まず、当選回数と当落の関係を図示する。

sp1 <- ggplot(HR09LDP, aes(previous, wlsmd)) +
    labs(x = "当選回数", y = "当落")
sp1 + geom_point()

これで、図14.11 (p.272) と同じ図ができるが、これだとほとんどの点が重なっているので、点を少しずらして（“jitter” して）みよう。

sp1 + geom_jitter(position = position_jitter(height = 0.05))

これで少しは点の散らばり方が見えるようになった。

これと同様に、選挙費用と当落の関係を図示してみよう。

sp2 <- ggplot(HR09LDP, aes(expm, wlsmd)) + labs(x = "選挙費用（100万円）", y = "当落")
sp2 + geom_jitter(position = position_jitter(height = 0.05))

これで、図14.12 (p.273) と同じような図ができた。

また、相関係数は、

with(HR09LDP, cor(previous, wlsmd))

## [1] 0.3538019

with(HR09LDP, cor(expm, wlsmd))

## [1] 0.07971141

である。

14.3.3 ロジスティック回帰式を推定する

小選挙区での当落を応答変数、過去の当選回数と選挙費用を説明変数とするロジスティック回帰式を推定しよう。

fit.ldp09 <- glm(wlsmd ~ previous + expm, data = HR09LDP,
                 family = binomial(link = "logit"))
summary(fit.ldp09)

## 
## Call:
## glm(formula = wlsmd ~ previous + expm, family = binomial(link = "logit"), 
##     data = HR09LDP)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9991  -0.6848  -0.4752  -0.4158   2.1276  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.365327   0.444358  -5.323 1.02e-07
## previous     0.269494   0.049971   5.393 6.93e-08
## expm        -0.003712   0.034280  -0.108    0.914
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 306.11  on 289  degrees of freedom
## Residual deviance: 271.67  on 287  degrees of freedom
## AIC: 277.67
## 
## Number of Fisher Scoring iterations: 4

これで、推定結果（表14.9 [p.274] と同じ結果）が示された。

14.3.4 ロジスティック回帰式の精度を求める

次に、予測の正解率を求めてみよう。

predicted.ldp <- predict(fit.ldp09, type = "response")
error.rate.ldp <- mean((predicted.ldp > 0.5 & HR09LDP$wlsmd == 0) | 
                           predicted.ldp < 0.5 & HR09LDP$wlsmd ==1)
1 - error.rate.ldp

## [1] 0.7655172

表14.10と同じく、予測の正解率は76.55%であることがわかった。

14.3.5 回帰係数を包括的に検定する

検定のために、\(\chi^2\)値と\(p\)値を求めよう。

(chisq.ldp <- with(fit.ldp09, null.deviance - deviance))

## [1] 34.4399

(df.ldp <- with(fit.ldp09, df.null - df.residual))

## [1] 2

pchisq(chisq.ldp, df.ldp, lower.tail = FALSE)

## [1] 3.322538e-08

これで、教科書 pp.274-275 と同じ結果が得られた。

14.3.6 回帰係数を個別に検定する

個別の検定は、係数の\(p\)値を使って行うことができる。すでに回帰式の推定は終わっているので、summary() で推定結果を表示すればよい。

summary(fit.ldp09)

## 
## Call:
## glm(formula = wlsmd ~ previous + expm, family = binomial(link = "logit"), 
##     data = HR09LDP)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9991  -0.6848  -0.4752  -0.4158   2.1276  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.365327   0.444358  -5.323 1.02e-07
## previous     0.269494   0.049971   5.393 6.93e-08
## expm        -0.003712   0.034280  -0.108    0.914
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 306.11  on 289  degrees of freedom
## Residual deviance: 271.67  on 287  degrees of freedom
## AIC: 277.67
## 
## Number of Fisher Scoring iterations: 4

この結果のうち“Pr(> |z|)” の列に示されている数字を読み、教科書と同じように検定を行えばよい。

14.3.7 説明変数を複数の値に固定したグラフを作成する

まず、表14.12 (p.277) を再現しよう。

prev.mean <- mean(HR09LDP$previous)  # previousの平均値
prev.sd <- sd(HR09LDP$previous)      # previousのsd
# previousの8つの値
prev.vec <- c(min(HR09LDP$previous), max(HR09LDP$previous), 0, 1,
              prev.mean - 1/2, prev.mean + 1/2,
              prev.mean - prev.sd/2, prev.mean + prev.sd/2)
expm.mean <- mean(HR09LDP$expm)  # expmnの平均値
expm.sd <- sd(HR09LDP$expm)      # expmのsd
## expmの8つの値
expm.vec <- c(min(HR09LDP$expm), max(HR09LDP$expm), 0, 1,
              expm.mean - 1/2, expm.mean + 1/2,
              expm.mean - expm.sd/2, expm.mean + expm.sd/2)
k <- 8  #　調べる値の数
wl.prev <- predict(fit.ldp09, type = "response",
                   newdata = list( previous = prev.vec, expm = rep(expm.mean, k)))
## expmの値を変えて当選確率を予測する：previousは平均値に固定
wl.expm <- predict(fit.ldp09, type = "response",
                   newdata = list(previous = rep(prev.mean, k), expm = expm.vec))
## 結果をまとめる
res <- rbind(wl.prev, wl.expm)
tbl.14.12 <- cbind(res[,2] - res[,1],
                  res[,4] - res[,3],
                  res[,6] - res[,5],
                  res[,8] - res[,7])
colnames(tbl.14.12) <- c("min -> max", "0 -> 1", "-+1/2", "-+sd/2")
tbl.14.12

##          min -> max        0 -> 1         -+1/2       -+sd/2
## wl.prev  0.78778462  0.0228538753  0.0420956804  0.129962524
## wl.expm -0.01464802 -0.0005940443 -0.0005796831 -0.002629427

これで、説明変数の変化に応じて応答変数がどのように変化するかがわかった。

候補者の選挙費用を最小値に固定した分析

選挙費用 (expm) が観測された最小値のとき、当選回数 (previous) の値の変化によって当選確率がどのように変化するかを示す図を作る。

まず、当選回数の変化に応じた予測値を求める。求める範囲は、当選回数 (previous) の観測された最小値と最大値の間とする。

HR09LDP <- mutate(HR09LDP,
                  prev.vec = seq(min(previous), max(previous), length = n()),
                  p.hat.min = predict(fit.ldp09, type = "response",
                                      list(previous = prev.vec, 
                                           expm = rep(min(expm), n()))))

これを図示してみよう。

wl.prev.lab <- labs(x = "当選回数", y = "当選確率")
p.hat.min <- ggplot(HR09LDP, aes(prev.vec, p.hat.min)) + geom_line(color = "red")
p.hat.min + wl.prev.lab + ggtitle("選挙費用が最小の場合") + ylim(0, 1)

これで、図14.13 と同じ図ができた。

候補者の選挙費用を最大値に固定した分析

選挙費用 (expm) が観測された最大値のも、同様に図示してみよう。

まず、当選回数の変化に応じた予測値を求める。求める範囲は、当選回数 (previous) の観測された最小値と最大値の間とする。

exp.max <- with(HR09LDP, max(expm))
HR09LDP <- mutate(HR09LDP,
                  p.hat.max = predict(fit.ldp09, type = "response",
                                      list(previous = prev.vec,
                                           expm = rep(max(expm), n()))))

これを図示してみよう。

p.hat.max <- ggplot(HR09LDP, aes(prev.vec, p.hat.max)) + geom_line(color = "blue")
p.hat.max + wl.prev.lab + ggtitle("選挙費用が最大の場合") + ylim(0, 1)

これで、図14.14 (p.281) と同じ図ができた。

最後に、2つの図を重ねて示そう。 ggplot2 で作図するための準備として、新たにP.hat.prev というデータフレームを作る。

P.hat.prev <- melt(HR09LDP, id = "prev.vec",
                   measure = c("p.hat.min", "p.hat.max"))
head(P.hat.prev)

##     prev.vec  variable      value
## 1 0.00000000 p.hat.min 0.08585220
## 2 0.05536332 p.hat.min 0.08703041
## 3 0.11072664 p.hat.min 0.08822323
## 4 0.16608997 p.hat.min 0.08943079
## 5 0.22145329 p.hat.min 0.09065325
## 6 0.27681661 p.hat.min 0.09189072

tail(P.hat.prev)

##     prev.vec  variable     value
## 575 15.72318 p.hat.max 0.8554332
## 576 15.77855 p.hat.max 0.8572686
## 577 15.83391 p.hat.max 0.8590844
## 578 15.88927 p.hat.max 0.8608810
## 579 15.94464 p.hat.max 0.8626583
## 580 16.00000 p.hat.max 0.8644165

このデータフレームを利用して作図しよう。

p.hat.ldp <- ggplot(P.hat.prev, aes(prev.vec, value)) + wl.prev.lab
p.hat.ldp <- p.hat.ldp + geom_line(aes(color = variable)) + ylim(0, 1)
p.hat.ldp + 
    scale_color_discrete(name="選挙費用（100万円）",
                         labels=c(paste0("最小値 (", min(HR09LDP$expm), ")"), 
                                  paste0("最大値 (", max(HR09LDP$expm), ")")))

これで、図14.15 (p.281) と同じ図ができた。

『Stataによる計量政治学』（浅野正彦, 矢内勇生)：Rによる分析

Yuki Yanai

October 16, 2014 (rev. June 6, 2015)

第14章　ロジスティック回帰分析

14.0 準備

14.1 ロジスティック関数

14.2 ロジスティック回帰分析の手順

14.2.1 帰無仮説と対立仮説を設定する

14.2.2 説明変数と応答変数の散布図

14.2.3 ロジスティック回帰式の推定値を求める

14.2.4 ロジスティック回帰式の精度を求める

14.2.5 回帰係数を包括的に検定する

14.2.6 回帰係数を個別に検定する

14.2.7 説明変数を複数の値に固定したグラフを作成する

当選0回の候補者に関する条件付き効果プロット

当選5回の候補者に関する条件付き効果プロット

14.3 衆院選データを使ったロジスティック回帰分析

14.3.1 帰無仮説と対立仮説

14.3.2 説明変数と応答変数の散布図

14.3.3 ロジスティック回帰式を推定する

14.3.4 ロジスティック回帰式の精度を求める

14.3.5 回帰係数を包括的に検定する

14.3.6 回帰係数を個別に検定する

14.3.7 説明変数を複数の値に固定したグラフを作成する

候補者の選挙費用を最小値に固定した分析

候補者の選挙費用を最大値に固定した分析

『Stataによる計量政治学』（浅野正彦, 矢内勇生)：Rによる分析

Yuki Yanai

October 16, 2014 (rev. June 6, 2015)

第14章 ロジスティック回帰分析

14.0 準備

14.1 ロジスティック関数

14.2 ロジスティック回帰分析の手順

14.2.1 帰無仮説と対立仮説を設定する

14.2.2 説明変数と応答変数の散布図

14.2.3 ロジスティック回帰式の推定値を求める

14.2.4 ロジスティック回帰式の精度を求める

14.2.5 回帰係数を包括的に検定する

14.2.6 回帰係数を個別に検定する

14.2.7 説明変数を複数の値に固定したグラフを作成する

当選0回の候補者に関する条件付き効果プロット

当選5回の候補者に関する条件付き効果プロット

14.3 衆院選データを使ったロジスティック回帰分析

14.3.1 帰無仮説と対立仮説

14.3.2 説明変数と応答変数の散布図

14.3.3 ロジスティック回帰式を推定する

14.3.4 ロジスティック回帰式の精度を求める

14.3.5 回帰係数を包括的に検定する

14.3.6 回帰係数を個別に検定する

14.3.7 説明変数を複数の値に固定したグラフを作成する

候補者の選挙費用を最小値に固定した分析

候補者の選挙費用を最大値に固定した分析

第14章　ロジスティック回帰分析