必要なパッケージを読み込む。
::p_load(tidyverse,
pacman broom)
次に、ggplot2のテーマとフォントの設定を行う。自分好みの設定がある場合は自由に変えてよい。
if (.Platform$OS.type == "windows") {
if (require(fontregisterer)) {
<- "Yu Gothic"
my_font else {
} <- "Japan1"
my_font
}else if (capabilities("aqua")) {
} <- "HiraginoSans-W3"
my_font else {
} <- "IPAexGothic"
my_font
}
theme_set(theme_gray(base_size = 9,
base_family = my_font))
説明のために『Rによる計量政治学』(浅野正彦, 矢内勇生. 2018)で使用されているデータ(hr-data.csv) を使う。 まず、このデータをダウンロードする。
download.file(
url = "https://raw.githubusercontent.com/yukiyanai/quant-methods-R/master/data_fixed/hr-data.csv",
dest = "data/hr-data.csv"
)
ダウンロードしたデータを読み込む。
<- read_csv("data/hr-data.csv") HR
小選挙区での当落を表すダミー変数(これを二値の応答変数として扱う)を作る。
<- HR %>%
HR mutate(wlsmd = ifelse(wl == "当選", 1, 0))
2009年の結果だけ抜き出し、HR09として保存する。
<- HR %>%
HR09 filter(year == 2009)
expが欠測の個体が1つあるので、(とりあえず)除外しておく。
<- na.omit(HR09) HR09
例として、2009年の衆院選データを用い、選挙費用(expm
, 単位100万円)で衆院選小選挙区での当落 (wlsmd
) を説明するモデルを考える。このモデルを式で表すと、 \[
\mathrm{Pr}(\text{当選}) = \mathrm{logit}^{-1}(\alpha + \beta \text{選挙費用})
\] となる。
まず、2つの変数の関係を図示してみよう。
<- ggplot(HR09, aes(x = expm, y = wlsmd)) +
p geom_jitter(size = 1,
alpha = 1/3,
width = 0,
height = 0.05) +
labs(x = "選挙費用(100万円)",
y = "小選挙区での当落")
plot(p)
これに(無理やり)通常の回帰直線を当てはめると、
plot(p + geom_smooth(method = "lm"))
となる。
Rでロジスティック回帰分析を行うにはglm()
を使う。 この関数は、一般化線形モデル (Generalized Linear Models: GLM) を当てはめるための関数で、ロジスティック回帰以外でも頻繁に使う関数である。 この関数をロジスティック回帰分析に使うには、引数family = binomial(link = "logit")
を指定する。
<- glm(wlsmd ~ expm,
fit1 data = HR09,
family = binomial(link = "logit"))
推定結果を確認する。
tidy(fit1, conf.int = TRUE)
## # A tibble: 2 × 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -2.01 0.126 -16.0 2.82e-57 -2.26 -1.77
## 2 expm 0.146 0.0142 10.3 1.14e-24 0.118 0.174
これで、推定はできた。
問題:推定値はどう解釈する?
次に、成功確率の予測値を求めよう。 予測値を求めるには、OLSのときと同様、predict()
を使えばよい。 ただし、predict()
のデフォルトでは、線形予測子の測定単位で予測値を返す。 ロジスティック回帰の線形予測子\(X\beta\) は対数オッズスケール \((-\infty, \infty)\) なので、 predict(fit1)
とすると、成功確率が返ってこない(各自で確かめてみること)。 私たちが欲しいのは確率なので、type = "response"
を指定し、応答変数の測定単位で値を得る(ここで、応答 [response] と考えられているのは、当落そのものではなく、当選確率であることに注意)。
<- HR09 %>%
HR09 mutate(pred = predict(fit1, type = "response"))
これを図示してみよう。
<- p +
p_pred geom_line(data = HR09,
aes(x = expm, y = pred),
color = "dodgerblue") +
labs(x = "選挙費用(100万円)",
y = "小選挙区での当選確率")
plot(p_pred)
あるいは、ggplot2の機能を利用して、次のようにすることもできる。
<- ggplot(HR09, aes(x = expm, y = wlsmd)) +
p3 geom_jitter(size = 1,
alpha = 1/3,
width = 0,
height = 0.05) +
geom_smooth(method = "glm",
method.args = list(family = binomial(link = "logit"))) +
labs(x = "選挙費用(100万円)",
y = "小選挙区での当選確率")
print(p3)
説明変数が増えても、基本的な分析法は変わらない。 先ほどまでのモデルに、過去の当選回数 (previous
) という説明変数を加えてみよう。
ロジスティック回帰を当てはめると、次のような結果が得られる。
<- glm(wlsmd ~ expm + previous,
fit2 data = HR09,
family = binomial(link = "logit"))
tidy(fit2, conf.int = TRUE)
## # A tibble: 3 × 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -2.03 0.132 -15.4 2.31e-53 -2.30 -1.78
## 2 expm 0.0500 0.0181 2.77 5.62e- 3 0.0144 0.0853
## 3 previous 0.324 0.0370 8.75 2.06e-18 0.253 0.399
問題:この推定結果をどう解釈する?
この結果を図示したいが、説明変数が複数あるので、先ほどのように簡単に図示することはできない。 選挙費用の効果は統計的に有意でない(ただし、符号は期待通り正)ので、選挙費用を平均値に固定し、過去の当選回数と当選確率の関係を図示してみよう。
選挙費用が平均値のときの当選確率の予測値は次のように求める。
<- HR09 %>%
HR09 mutate(pred2 = predict(fit2,
type = "response",
list(previous = previous,
expm = rep(mean(expm), n()))))
これを図示しよう。
<- ggplot(HR09, aes(x = previous)) +
p_pred2 geom_jitter(aes(y = wlsmd),
size = 1,
alpha = 1/3,
width = 0.1,
height = 0.05) +
geom_line(aes(y = pred2), color = "tomato") +
labs(x = "当選回数",
y = "小選挙区での当選確率")
print(p_pred2)
これに信頼区間を加えよう。そのために、まず各予測値の標準誤差を求める。
<- HR09 %>%
HR09 mutate(se2 = predict(fit2, type = "response",
list(previous = previous,
expm = rep(mean(expm), n())),
se.fit = TRUE)$se.fit)
予測値$$2seの範囲を図示する(つまり、ほぼ95%信頼区間を描く)ことにして、その範囲を求める。
<- HR09 %>%
HR09 mutate(fit2_lb = pred2 - 2 * se2,
fit2_ub = pred2 + 2 * se2)
求めた区間を図示しよう。
<- p_pred2 +
p_pred2se geom_ribbon(data = HR09,
aes(ymin = fit2_lb, ymax = fit2_ub),
alpha = 0.2)
plot(p_pred2se)
ロジスティック回帰の結果に信頼区間を加えることができた。