Chapter 10 【発展】ロジスティック回帰分析（作成中）
本章では、前章で学んだ回帰分析を発展させて、複数の独立変数を扱う方法について説明する。
内容に入る前に、右上のプロジェクトのボックスの横が、前章で作成したプロジェクトの名前（たとえば、seminar_sociology_r）になっているかどうかを確認しておこう。なっていない場合は、右上のボックスをクリックして、「Open Project…」を選択し、前章で作成したRprojファイル（たとえば、seminar_sociology_r.Rprojといったような名前になっている）を選んで、プロジェクトを切り替えよう。
さらに、これまでの章で説明した以下のパッケージを読み込んだ上で、第4章で作成したデータを読み込んでpiaacというデータフレームに入れていることを前提とする。具体的には、以下のコードを実行しておく必要がある。
library(tidyverse)
library(gtsummary)
library(flextable)
library(modelsummary)
library(ggeffects) #第9章（交互作用）で紹介
piaac <- read_rds("data/piaac_sample_analytic.rds")第5章で確認したように、ggplotの設定を変更しておくことで見やすいグラフを作ることができる。ここでは以下のコードを実行している。
Macの場合：
theme_set(theme_bw(
base_family = "HiraginoSans-W3",
base_size = 11,
base_rect_size = 0.2,
base_line_size = 0.2
))Windowsの場合：
10.1 二値変数を従属変数にする：線形回帰分析を使ったアプローチ
10.1.1 線形確率モデル
従属変数が2値のカテゴリ変数（0/1）の場合ですでにみたように、0/1の値をとる二値変数を従属変数として回帰分析を使うことができる。
\[ y = \beta_0 + \beta_1x_1 + \cdots + \beta_kx_k \]
このとき、係数\(\beta_1\)は他の変数\(x_2, \cdots, x_k\)を一定としたうえで、\(x_1\)が1単位高いと\(y\)（従属変数が1をとる割合）がどれだけ高いのかを表す。
たとえば、この1年間に職場での訓練（OJTとよぶ）を受ける率が性別によってどの程度異なるのかを知りたいとしたら、次のような式を推定する。
##
## Call:
## lm(formula = ojt ~ gender, data = piaac)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.4079 -0.4079 -0.3263 0.5921 0.6737
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.32632 0.01343 24.30 < 2e-16 ***
## gender男性 0.08157 0.01846 4.42 1.03e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4811 on 2726 degrees of freedom
## Multiple R-squared: 0.007116, Adjusted R-squared: 0.006751
## F-statistic: 19.54 on 1 and 2726 DF, p-value: 1.026e-05また、年齢とOJT受講の関係が線形であると仮定して、年齢および学歴が同じであったとしてもなお、性別によってOJTを受ける割合が異なっているかどうかを知りたいとしたら、次のような式を推定する。
##
## Call:
## lm(formula = ojt ~ gender + age + educ, data = piaac)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5679 -0.3590 -0.2487 0.4981 0.8575
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.2790638 0.0539055 5.177 2.42e-07 ***
## gender男性 0.0532179 0.0188275 2.827 0.00474 **
## age -0.0021335 0.0008478 -2.516 0.01192 *
## educ高校 0.0635070 0.0349841 1.815 0.06959 .
## educ短大高専 0.1652651 0.0369943 4.467 8.24e-06 ***
## educ大学大学院 0.2890064 0.0356772 8.101 8.17e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.469 on 2722 degrees of freedom
## Multiple R-squared: 0.05794, Adjusted R-squared: 0.05621
## F-statistic: 33.48 on 5 and 2722 DF, p-value: < 2.2e-16この結果から、年齢および学歴が同程度であったとしても、男性がOJTを受ける割合は女性と比べて0.053ポイント（5.3 %ポイント）高いといえる。
このように、0/1の値をとる二値変数を従属変数として回帰分析を推定するような場合を指して、線形確率モデル（linear probability model）と呼ぶ。
10.1.2 残差が正規分布しない問題への対処
回帰分析で係数の標準誤差を求めるときには、残差が正規分布するという仮定が置かれている。しかしながら、線形確率モデルの場合はこの仮定は通常満たされない。このことを、不均一分散 heteroskedasticityという。不均一分散が生じている場合には、標準誤差にバイアスが生じる。そこで、不均一分散に対して頑健な標準誤差（ロバスト標準誤差 robust standard error）を用いることで、適切な標準誤差を計算することができる。
ロバスト標準誤差を求める場合には、estimatrパッケージに含まれているlm_robust()関数を使うのが手っ取り早い。通常のlm()を、lm_robust()に変えるだけで、ロバスト標準誤差を計算してくれる。
lm()で求めた推定結果と、lm_robust()で求めた推定結果を比べてみよう：
| 通常のSE | Robust SE | |
|---|---|---|
| (Intercept) | 0.2791*** | 0.2791*** |
| (0.0539) | (0.0518) | |
| gender男性 | 0.0532** | 0.0532** |
| (0.0188) | (0.0188) | |
| age | -0.0021* | -0.0021* |
| (0.0008) | (0.0008) | |
| educ高校 | 0.0635+ | 0.0635* |
| (0.0350) | (0.0308) | |
| educ短大高専 | 0.1653*** | 0.1653*** |
| (0.0370) | (0.0338) | |
| educ大学大学院 | 0.2890*** | 0.2890*** |
| (0.0357) | (0.0326) | |
| Num.Obs. | 2728 | 2728 |
| R2 | 0.058 | 0.058 |
| R2 Adj. | 0.056 | 0.056 |
| AIC | 3618.6 | 3618.6 |
| BIC | 3660.0 | 3660.0 |
| Log.Lik. | -1802.312 | |
| F | 33.484 | |
| RMSE | 0.47 | 0.47 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | ||
通常の標準誤差を用いた場合と、ロバスト標準誤差を用いた場合とでは、わずかに標準誤差の値が異なっていることがわかる。このずれがどれくらい大きくなるかは、データや変数の性質によって変わってくる。ただし、係数の値には変化はないため、実質的な結論に違いが出ることはあまりない。線形確率モデルを使う場合には、基本的にはロバスト標準誤差を使うことが推奨されている。
10.2 ロジスティック回帰分析とは
10.2.1 ロジスティック回帰分析（ロジットモデル）
\[ \log \frac{\Pr(y = 1)}{1 - \Pr(y = 1)} = \beta_0 + \beta_1x_1 + \cdots + \beta_kx_k \] または、式を変形して次のように書く。
\[ \begin{align} \log \frac{\Pr(y = 1)}{1 - \Pr(y = 1)} &= \beta_0 + \beta_1x_1 + \cdots + \beta_kx_k \\ \frac{\Pr(y = 1)}{1 - \Pr(y = 1)} &= \exp(\beta_0 + \beta_1x_1 + \cdots + \beta_kx_k) \\ \Pr(y = 1) &= (1 - \Pr(y = 1))\exp(\beta_0 + \beta_1x_1 + \cdots + \beta_kx_k) \\ (1 + \exp(\beta_0 + \beta_1x_1 + \cdots + \beta_kx_k))\Pr(y = 1) &= \exp(\beta_0 + \beta_1x_1 + \cdots + \beta_kx_k) \\ \Pr(y = 1) &= \frac{\exp(\beta_0 + \beta_1x_1 + \cdots + \beta_kx_k)}{1 + \exp(\beta_0 + \beta_1x_1 + \cdots + \beta_kx_k)} \end{align} \]
##
## Call:
## glm(formula = ojt ~ gender, family = binomial, data = piaac)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.72486 0.05952 -12.178 < 2e-16 ***
## gender男性 0.35218 0.08006 4.399 1.09e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 3593.8 on 2727 degrees of freedom
## Residual deviance: 3574.3 on 2726 degrees of freedom
## AIC: 3578.3
##
## Number of Fisher Scoring iterations: 4「gender男性」の行は、女性と比べて男性であると対数オッズ（log-odds）がcoef高いことを表している。
10.2.2 対数オッズ比の解釈
急に出てきた「対数オッズ比」とはなんだろうか？これを理解するために、クロス集計表を作ってみよう。性別ごとにOJTを受ける・受けない人数を集計したクロス集計表は次のようになる。ここでは説明を分かりやすくするために順序を変更している。
piaac %>%
mutate(gender = fct_relevel(gender, "男性", "女性")) %>%
mutate(ojt = factor(ojt, levels = 1:0, labels = c("1（受けた）", "0（受けなかった）"))) %>%
tbl_cross(gender, ojt)| ojt | Total | ||
|---|---|---|---|
| 1（受けた） | 0（受けなかった） | ||
| gender | |||
| 男性 | 589 | 855 | 1,444 |
| 女性 | 419 | 865 | 1,284 |
| Total | 1,008 | 1,720 | 2,728 |
この表から、「OJTを受けなかった人数に対する受けた人数の比」を計算してみる。
男性の場合はこのように：
\[ 589 / 855 = 0.689 \]
女性の場合はこのように：
\[ 419 / 865 = 0.484 \]
それぞれ計算できる。このように、0をとる人数（確率）に対する1をとる人数（確率）の比を取った値のことを指して、オッズ（odds）という。オッズは1をとる人数と0をとる人数がちょうど同じときに1となり、1をとる人数（分子）が0をとる人数（分母）に対して少ないほど0に近づいていき、逆に1をとる人数（分子）が0をとる人数（分母）に対して多いほど無限大に近づいていく。
つぎに、こうして計算したオッズどうしの比（オッズ比 odds ratio）をとった値を計算してみよう。
\[ \frac{589 / 855}{419 / 865} = 1.422 \]
この値は、「女性のオッズに対して男性のオッズが何倍であるか」を示しており、1をとるときには女性と男性のオッズが同じであることを意味し、1よりも大きいときには女性のオッズ（分母）よりも男性のオッズ（分子）のほうが大きいことを意味し、1よりも小さいときには男性のオッズ（分子）よりも女性のオッズ（分母）のほうが大きいことを意味する。
このオッズ比の自然対数をとったのが、対数オッズ比 log-odds ratioである。
\[ \log \left( \frac{589 / 855}{419 / 865} \right) = 0.352 \]
この値は、「女性の対数オッズに対して男性の対数オッズがどれくらい大きいか」を示しており、0をとるときには女性と男性の（対数）オッズが同じであることを意味し、0よりも大きいときには女性の（対数）オッズ（分母）よりも男性の（対数）オッズ（分子）のほうが大きいことを意味し、0よりも小さいときには男性の（対数）オッズ（分子）よりも女性の（対数）オッズ（分母）のほうが大きいことを意味する。
この対数オッズ比の値が、先ほど推定したロジスティック回帰分析の「gender男性」の係数の値に一致していることを確認しよう。このことが、ロジスティック回帰分析の係数が対数オッズ比を表すという意味である。
今までの話を一般化しよう。次のような2×2のクロス集計表があるとする。
| Y = 1 | Y = 0 | |
|---|---|---|
| X = 1 | \(N_{11}\) | \(N_{10}\) |
| X = 0 | \(N_{01}\) | \(N_{00}\) |
このとき、オッズ比および対数オッズ比はそれぞれ次のように定義され、次のような特徴をもつ。
| オッズ比 \((N_{11} / N_{10}) / (N_{01} / N_{00})\) | 対数オッズ比 \(\log((N_{11} / N_{10}) / (N_{01} / N_{00}))\) | |
|---|---|---|
| X = 0と比べてX = 1のほうがよりY = 1をとりやすい場合 | 1よりも大きくなり、大きくなるほど∞に近づく | 0よりも大きくなり、大きくなるほど+∞に近づく |
| Y = 1のとりやすさが等しい場合 | 1をとる | 0をとる |
| X = 1と比べてX = 0のほうがよりY = 1をとりやすい場合 | 1よりも小さくなり、小さくなるほど0に近づく | 0よりも小さくなり、小さくなるほど-∞に近づく |
なお、対数の性質より、対数オッズ比は\(\log((N_{11} / N_{10}) / (N_{01} / N_{00})) = \log(N_{11} / N_{10}) - (N_{01} / N_{00})\)などと書くこともできる。
10.2.3 連続変数を独立変数とする場合
連続変数を独立変数とする場合は先のような表をイメージするのは難しくなってくるが、解釈は類似している。たとえば次のように性別と年齢を独立変数とするロジスティック回帰分析を推定する。
##
## Call:
## glm(formula = ojt ~ gender + age, family = binomial, data = piaac)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.006992 0.171320 -0.041 0.967
## gender男性 0.350124 0.080351 4.357 1.32e-05 ***
## age -0.016442 0.003705 -4.437 9.10e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 3593.8 on 2727 degrees of freedom
## Residual deviance: 3554.5 on 2725 degrees of freedom
## AIC: 3560.5
##
## Number of Fisher Scoring iterations: 4このときの年齢の係数は年齢が1高いとOJTを受ける対数オッズが-0.016だけ高いということを意味している。
10.3 平均限界効果
ロジスティック回帰分析では係数は対数オッズを表す。しかし、対数オッズ比は線形確率モデルと比べると直感的に解釈しにくい。そこで、ロジスティック回帰分析を使う場合には、平均限界効果 Average marginal effect, AMEと合わせて結果を解釈することが強く推奨される。
平均限界効果とはなんだろうか。これを理解するため、先ほどのロジスティック回帰分析を次のように変形する（式変形については対数の知識が必要）：
## # A tibble: 2,728 × 3
## id gender age
## <int> <fct> <dbl>
## 1 1 女性 45
## 2 2 男性 48
## 3 3 男性 47
## 4 4 男性 58
## 5 6 女性 35
## 6 9 男性 56
## 7 12 女性 32
## 8 15 男性 45
## 9 18 男性 30
## 10 19 男性 46
## # ℹ 2,718 more rowsreg_res3 <- lm(data = piaac, ojt ~ gender + age + educ)
logit_res3 <- glm(data = piaac, ojt ~ gender + age + educ, family = binomial())
margins_logit_res3 <- margins(logit_res3)modelsummary(list("LPM" = reg_res3, "Logit" = logit_res3, "AME" = margins_logit_res3),
stars = TRUE)| LPM | Logit | AME | |
|---|---|---|---|
| (Intercept) | 0.279*** | -1.022*** | |
| (0.054) | (0.258) | ||
| gender男性 | 0.053** | 0.247** | 0.054** |
| (0.019) | (0.086) | (0.019) | |
| age | -0.002* | -0.010* | -0.002* |
| (0.001) | (0.004) | (0.001) | |
| educ高校 | 0.064+ | 0.360* | 0.066* |
| (0.035) | (0.182) | (0.031) | |
| educ短大高専 | 0.165*** | 0.832*** | 0.169*** |
| (0.037) | (0.188) | (0.034) | |
| educ大学大学院 | 0.289*** | 1.329*** | 0.290*** |
| (0.036) | (0.181) | (0.033) | |
| Num.Obs. | 2728 | 2728 | |
| R2 | 0.058 | ||
| R2 Adj. | 0.056 | ||
| AIC | 3618.6 | 3445.4 | 3445.4 |
| BIC | 3660.0 | 3480.8 | 3480.8 |
| Log.Lik. | -1802.312 | -1716.687 | |
| F | 33.484 | 30.326 | |
| RMSE | 0.47 | 0.47 | 0.47 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | |||