第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 devianceResidual 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 のそれぞれについて、調べる対象となる値のベクトルを作る。 必要な値は、

  1. 最小値
  2. 最大値
  3. 0
  4. 1
  5. 平均値 - 1/2
  6. 平均値 + 1/2
  7. 平均値 - sd/2
  8. 平均値 + 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))