第12回講義資料

ロジスティック回帰分析

スライド

新しいタブで開く

セットアップ

 本日の実習用データ(LMSからダウンロード可能)と必要なパッケージ({tidyverse}、{summarytools}、{marginaleffects})を読み込む。データを読み込み、dfという名のオブジェクトとして格納する。

Code 01
library(tidyverse)
library(summarytools)
library(marginaleffects)

df <- read_csv("Data/Hino_Song.csv")

df

 各変数の詳細はスライドを参照すること。データ分析を進める前に、まず{summarytools}のdescr()関数を使ってdfの各変数の記述統計量を確認する。今回のデータの場合、全て連続変数(ダミー変数も連続変数として扱う)であるため、dfをそのままdescr()に渡せば良い。

Code 02
df %>%
  descr(stats = c("mean", "sd", "min", "max", "n.valid"),
        transpose = TRUE, order = "p")
Descriptive Statistics  
df  
N: 2555  

                   Mean   Std.Dev     Min     Max   N.Valid
--------------- ------- --------- ------- ------- ---------
          Voted    0.71      0.45    0.00    1.00   2555.00
         Female    0.46      0.50    0.00    1.00   2555.00
            Age   50.91     15.69   18.00   86.00   2555.00
           Educ    3.21      0.90    1.00    4.00   2555.00
      Knowledge    3.95      1.12    1.00    5.00   2555.00
       Ideology    5.40      2.16    0.00   10.00   2555.00

 今回は以下の問いに答えるモデルの推定を行う。

 有権者の投票参加を規定する要因を調べたい。投票所に足を運ぶには予め投票先を決めておく必要があろう。しかし、数多い選択肢(候補者、政党)の中から自分の望みを実現してくれそうな選択肢を見つけることは簡単な作業ではない。政治に関する知識があれば、比較的簡単に見つかるため、投票参加しやすいと考えられる。一方、そうでない有権者は自分にとっても最適な選択肢を見つけることを諦め、棄権するだろう。これは本当だろうか。

 まず、応答変数は回答者の投票参加の有無(Voted)である。こちらは投票すれば1、危険すれば0の値を取るダミー変数であり、応答変数として使う場合は二値変数、またはバイナリー変数とも呼ばれる。続いて、主な説明変数は回答者の主観的な政治知識(Knowledge)である。また、主な説明変数以外の説明変数(=統制変数)として性別(Female)、年齢(Age)、学歴(Educ)、イデオロギー(Ideology)を投入する。このモデルを可視化すると以下のようになる。

線形確率モデル

 まずは、これまで紹介してきた線形回帰分析を使って、モデルを推定してみよう。推定する回帰式は以下の通りである。

\[ \widehat{\mbox{Voted}} = \alpha + \beta_1 \mbox{Knowledge} + \beta_2 \mbox{Female} + \beta_3 \mbox{Age} + \beta_4 \mbox{Educ} + \beta_5 \mbox{Ideology} \]

Code 03
lm_fit <- lm(Voted ~ Knowledge + Female + Age + Educ + Ideology, data = df)

summary(lm_fit)

Call:
lm(formula = Voted ~ Knowledge + Female + Age + Educ + Ideology, 
    data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.0057 -0.4242  0.1452  0.2844  0.8522 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.1093644  0.0558137  -1.959 0.050168 .  
Knowledge    0.1186971  0.0077281  15.359  < 2e-16 ***
Female      -0.0549171  0.0170946  -3.213 0.001332 ** 
Age          0.0046902  0.0005579   8.408  < 2e-16 ***
Educ         0.0357970  0.0096391   3.714 0.000209 ***
Ideology     0.0040744  0.0038653   1.054 0.291944    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.4202 on 2549 degrees of freedom
Multiple R-squared:  0.1451,    Adjusted R-squared:  0.1435 
F-statistic: 86.55 on 5 and 2549 DF,  p-value: < 2.2e-16

 推定の結果、政治知識(Knowledge)の係数は0.119が得られた。これは「政治知識が1単位上がると、投票参加の確率は約11.9%ポイント上がる」ことを意味する(「%」でなく、「%ポイント」であることに注意)。\(p\)値も非常に小さく(\(p < 0.05\))、統計的に有意な関係であると言えよう。

 このようにバイナリー変数を応答変数として使う線形回帰モデルを線形確率モデル(linear probability model; LPM)と呼ぶ。これは非常に直感的な方法であるものの、一つ大きな問題がある。たとえば、主観的政治知識(Knowledge = 5)が5、男性(Female = 0)、86歳(Age = 86)、教育水準が4(Educ = 4)、イデオロギーが10(Ideology = 10)の回答者がいる場合、投票に参加する確率の予測値を計算してみよう。

Code 04
-0.1093644 + 0.1186971 * 5 - 0.0549171 * 0 + 0.0046902 * 86 + 0.0357970 * 4 + 0.0040744 * 10
[1] 1.07141

 主観的政治知識が5、男、86歳、教育水準が4、イデオロギーが10の回答者が投票に参加する確率は約107.1%である。確率は0〜100%の値を取るはずなのに、100%を超えてしまう。 他にも線形確率モデルは、線形回帰モデルの重要な仮定である分散の不均一性の満たさないなど、いくつかの問題がある。ここで登場するのがロジスティック回帰分析だ。

ロジスティック回帰分析

ロジスティック関数

 ロジスティック回帰分析を理解するためには、その背後にはるロジスティック関数に関する知識が必要である。ロジスティック関数は以下のような関数である。

\[ \text{logistic}(x) = \frac{1}{1 + e^{-x}} \]

 式内の\(e\)はネイピア数であり、2.71828184590…の無理数である。\(x\)\(-\infty \sim \infty\)の値を取り得る。たとえば、\(x\)\(-\infty\)の場合、\(e^{-x}\)\(e^{-(-\infty)} = \infty\)\(x\)\(\infty\)の場合、\(e^{-x}\)は$e^{-()} =\(0となる。つまり、\)e^{-x}$は\(0 \sim \infty\)の値を取る。この\(e^{-x}\)が0の場合、\(\text{logistic}(x)\)は1、\(e^{-x}\)\(\infty\)の場合、\(\text{logistic}(x)\)は0を取る。つまり、\(\text{logistic}(x)\)は0以上、1以下の取ることが分かる。確率もまた0以上1以下であるため、この関数が何らかの役割を果たすこととなる。以下の図は\(x\)の値ごとのロジスティック関数の折れ線グラフである。

 \(x\)の値が大きければ大きいほど、logistic(\(x\))の値は1へ近づくことが分かる。この\(x\)\(y^*\)と表記した場合、ロジスティック回帰分析は

\[ \begin{align} \mbox{Pr}(\mbox{Voted} = 1) & = \frac{1}{1 + e^{-y^*}} \\ y^* & = \alpha + \beta_1 \mbox{Knowledge} + \beta_2 \mbox{Female} + \beta_3 \mbox{Age} + \beta_4 \mbox{Educ} + \beta_5 \mbox{Ideology} \end{align} \]

における\(\alpha\)\(\beta_1\)、…を推定することである。ここでの\(y^*\)線形予測子(linear predictor)と呼ばれる。

実装

 ロジスティックの実装はglm()関数を使用し、使い方はlm()関数とほぼ同様である。

glm(応答変数 ~ 説明変数, data = データ名, family = binomial("logit"))

 lm()との違いはfamily引数が必要という点だ。ロジスティック回帰分析の場合はfamily = binomial("logit")、本講義では取り上げないものの同目的の分析手法であるプロビット回帰分析の場合はfamily = binomial("probit")を使用する。それでは以下のモデルを推定し、fit1に格納してみよう。

\[ \mbox{Pr}(\mbox{Voted} = 1) = \frac{1}{1 + e^{-(\alpha + \beta_1 \mbox{Knowledge} + \beta_2 \mbox{Female} + \beta_3 \mbox{Age} + \beta_4 \mbox{Educ} + \beta_5 \mbox{Ideology})}} \]

Code 05
fit1 <- glm(Voted ~ Knowledge + Female + Age + Educ + Ideology, 
            data = df, family = binomial("logit"))

summary(fit1)

Call:
glm(formula = Voted ~ Knowledge + Female + Age + Educ + Ideology, 
    family = binomial("logit"), data = df)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.3275  -0.9980   0.5583   0.7739   2.0312  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.30579    0.32049 -10.315  < 2e-16 ***
Knowledge    0.59290    0.04294  13.808  < 2e-16 ***
Female      -0.31597    0.09628  -3.282 0.001031 ** 
Age          0.02651    0.00320   8.284  < 2e-16 ***
Educ         0.20827    0.05430   3.835 0.000125 ***
Ideology     0.02392    0.02254   1.061 0.288586    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 3078.9  on 2554  degrees of freedom
Residual deviance: 2700.6  on 2549  degrees of freedom
AIC: 2712.6

Number of Fisher Scoring iterations: 4
ロジスティック回帰分析の決定係数?

ロジスティック回帰分析の場合、決定係数(\(R^2\))は表示されない。類似した概念として「疑似決定係数(Pseudo\(R^2\))」があるが、あまり使われない。どうしても疑似決定係数を出したい場合は{DescTools}パッケージのPseduoR2()関数を使う。

Code 06
# 予め{DescTools}をインストールしておくこと (コンソール上でinstall.packages("DescTools"))
DescTools::PseudoR2(fit1) # McFaddenの疑似決定係数
 McFadden 
0.1228484 

 それでは結果を解釈してみよう。線形確率モデルと同様、政治知識の\(p\)値は0.001未満であり、「主観的な政治知識と投票参加の間には統計的に有意な関係がある」ことが分かる。また、係数の符号(今回は+)を考えると、「主観的な政治知識が高くなると、投票に参加する確率も上がる」とも言えよう。しかし、政治知識が1単位上がると、投票参加の確率は具体的にどれくらい上がるか。線形確率モデルの場合は係数をそのまま解釈するだけで、具体的な上がりの度合いが分かったが、ロジスティック回帰分析では違う。

 たとえば、政治知識の係数は約0.593であるが、「主観的な政治知識が1上がると、投票参加の確率が0.593%p上がる」ということは完全な間違いである。これをどう解釈すべきだろうか。

 たおてば、主観的政治知識が3(Knowledge = 3)、女性(Female = 1)、20歳(Age = 20)、学歴が大卒(Educ = 4)、イデオロギーが中道(Ideology = 5)の場合の投票参加の予測確率を計算してみよう。まず、coef()関数で係数のみを抽出してみよう。

Code 07
fit1_coef <- coef(fit1) # coef()関数で係数のみを抽出する
fit1_coef
(Intercept)   Knowledge      Female         Age        Educ    Ideology 
-3.30579506  0.59290164 -0.31596873  0.02650943  0.20827008  0.02391612 

 続いて、切片には1を、その他の係数には各変数の具体的な値をかけて合計を求める。これが線形予測子だ。

Code 08
fit1_coef[1] * 1 + fit1_coef[2] * 3 + fit1_coef[3] * 1 + fit1_coef[4] * 20 + fit1_coef[5] * 4 + fit1_coef[6] * 5
(Intercept) 
 -0.3602094 

 同じ長さのベクトルが2つあれば、同じ位置同士の要素が計算されることを考えると、以下のように線形予測子を計算することもできる。

Code 09
sum(coef(fit1) * c(1, 3, 1, 20, 4, 5))
[1] -0.3602094

 あとは線形予測子の値をロジスティック関数に代入するだけだ。

Code 10
1 / (1 + exp(-(-0.3602094))) # 予測確率の計算
[1] 0.4109089

 結果は0.411であり、これは投票参加の予測確率は約41.1%であるこを意味する(以下の図を参照)。

 このような作業をより簡単にしてくれるのが、前回の講義で登場した{marginaleffects}パッケージである。このパッケージが提供するpredictions()関数の利用すれば予測確率が簡単に計算できる。使い方は以下の通りである。

predictions(回帰分析オブジェクト名, 
            newdata = datagrid(説明変数名1 = 値,
                               説明変数名2 = 値,
                               ...))

 たとえば、先程の例をpredictions()関数で計算してみよう。

Code 11
predictions(fit1, newdata = datagrid(Knowledge = 3, 
                                     Female    = 1,
                                     Age       = 20, 
                                     Educ      = 4, 
                                     Ideology  = 5))
  rowid     type predicted  std.error statistic     p.value  conf.low conf.high
1     1 response 0.4109089 0.02923061  14.05749 6.93037e-45 0.3550561 0.4691564
     Voted Knowledge Female Age Educ Ideology
1 0.709589         3      1  20    4        5

 datagrid()内で変数を特定の値で指定しない場合、その変数は自動的に平均値が代入される。たとえば、Knowledgeが1〜5の場合の投票参加確率の予測値を計算してみよう。その他の変数はすべて平均値に固定する。計算結果はfit1_predに格納する。

Code 12
fit1_pred <- predictions(fit1, newdata = datagrid(Knowledge = 1:5))

fit1_pred
  rowid     type predicted   std.error statistic       p.value  conf.low
1     1 response 0.3293996 0.027936256  11.79111  4.337629e-32 0.2771235
2     2 response 0.4705341 0.022008655  21.37950 2.073214e-101 0.4277233
3     3 response 0.6165432 0.013584204  45.38677  0.000000e+00 0.5895934
4     4 response 0.7441789 0.009492179  78.39916  0.000000e+00 0.7251353
5     5 response 0.8403326 0.009821140  85.56365  0.000000e+00 0.8201318
  conf.high    Voted   Female      Age     Educ Ideology Knowledge
1 0.3862676 0.709589 0.458317 50.91155 3.207436 5.400391         1
2 0.5137829 0.709589 0.458317 50.91155 3.207436 5.400391         2
3 0.6427954 0.709589 0.458317 50.91155 3.207436 5.400391         3
4 0.7623355 0.709589 0.458317 50.91155 3.207436 5.400391         4
5 0.8586556 0.709589 0.458317 50.91155 3.207436 5.400391         5

 datagrid()内で指定子なかったFemaleAgeEducIdeology列の値が最初に計算して記述統計量(の平均値)と一致することが分かる。今回はすべての列が必要なわけではないので、一部の列のみを出力してみよう。

Code 13
fit1_pred %>%
  select(Knowledge, predicted, std.error, conf.low, conf.high)
  Knowledge predicted   std.error  conf.low conf.high
1         1 0.3293996 0.027936256 0.2771235 0.3862676
2         2 0.4705341 0.022008655 0.4277233 0.5137829
3         3 0.6165432 0.013584204 0.5895934 0.6427954
4         4 0.7441789 0.009492179 0.7251353 0.7623355
5         5 0.8403326 0.009821140 0.8201318 0.8586556

 以上の結果を見れば、KnowledgeVotedに与える影響は一定ではないことがわかる。たとえば、Knowledgeの値が1の場合、投票に参加する確率は32.94%である。ここでKnowledgeの値が2の場合、投票に参加する確率は47.05%であり、14.11%p増加したことが分かる。それではKnowledgeの値が3の場合はどうだろうか。この場合の投票参加の確率は61.65%であり、先より14.60%p増加した。また、Knowledgeの値が4の場合の予測確率は74.42%(12.77%p増加)、5の場合のそれはは84.03%で、9.61%p増加したことが分かる。このようにロジスティック関数(から得られた曲線)の形が非線形であることを考えると、ある意味、当たり前のことであろう。

 したがって、ロジスティック回帰分析を行う場合、「Xが1単位上がるとYはZZ上がる/下がる」という表現はあまり使わず、文章では「正(負)の関係がある」、「統計的に有意な関係は見られない」程度とし、予測確率をグラフとして示すのが一般的である。たとえば、以上の結果を可視化したものが以下の図である。可視化については次回の講義で解説する。

Code 14
コード
fit1_pred %>%
  ggplot() +
  geom_pointrange(aes(x = Knowledge, y = predicted, 
                      ymin = conf.low, ymax = conf.high)) +
  labs(x = "主観的な政治知識", 
       y = "投票に参加する確率の予測値") +
  theme_bw(base_size = 14)

 Ageのように説明変数の中でも取り得る値が多い場合は以下のような図になってしまい、人によっては少し気持ち悪い図になってしまう。

Code 14
コード
predictions(fit1, newdata = datagrid(Age = 18:86)) %>%
  ggplot() +
  geom_pointrange(aes(x = Age, y = predicted, 
                      ymin = conf.low, ymax = conf.high)) +
  labs(x = "年齢 (歳)", y = "投票に参加する確率の予測値") +
  scale_x_continuous(breaks = c(20, 30, 40, 50, 60, 70, 80), 
                     labels = c(20, 30, 40, 50, 60, 70, 80)) +
  theme_bw(base_size = 14)

 したがって、取り得る値が多い連続変数を横軸にする場合は、折れ線グラフ (geom_line())とリボン (geom_ribbon())を併用することが一般的だ。

Code 14
コード
fit1_pred2 <- predictions(fit1, newdata = datagrid(Age = 18:86))

fit1_pred2 %>%
  ggplot() +
  geom_ribbon(aes(x = Age, ymin = conf.low, ymax = conf.high), 
              fill = "gray80") + 
  geom_line(aes(x = Age, y = predicted)) +
  labs(x = "年齢 (歳)", y = "投票に参加する確率の予測値") +
  scale_x_continuous(breaks = c(20, 30, 40, 50, 60, 70, 80), labels = c(20, 30, 40, 50, 60, 70, 80)) +
  theme_bw(base_size = 14)