回帰モデルの結果を限界効果（marginal effect）で示す

前回はロジスティック回帰モデルを（少しだけ）学び直してみました。

necostat.hatenablog.jp

ここで扱ったNorton先生の文献では、オッズ比に代わるロジスティック回帰モデルの指標として限界効果（marginal effect）を推奨していましたので、簡単にどんなものなのかとRでの実装方法をまとめてみます。

今回の参考文献は以下のものです。

Norton EC, Dowd BE, Maciejewski ML. Marginal Effects-Quantifying the Effect of Changes in Risk Factors in Logistic Regression Models. JAMA. 2019 Apr 2;321(13):1304-1305.

限界効果（marginal effect）とは
- 線形回帰モデルの限界効果
- ロジスティック回帰モデルの限界効果
限界効果の良いところ
- オッズ比は全範囲で発生確率を適切に反映しない
- 限界効果の方がモデルに依存しにくい
限界効果を使うときの注意点
- 想定する対象者を決めなくてはならない
- ベースラインの割合と比べる必要あり
おわりに
参考資料

限界効果（marginal effect）とは

説明変数が1単位増加するときに、アウトカムの期待値（2値アウトカムなら発生割合）がどれだけ変化するかを表したものです。説明変数が連続変数のときはmarginal effect、カテゴリー変数のときはincremental effectと呼ばれますが、まとめてmarginal effectと呼ばれることが多いようです。説明変数がカテゴリー変数の場合は比較的シンプルだと思うので、以後は説明変数が連続変数の場合の限界効果を考えます。

線形回帰モデルの限界効果

線形回帰モデルではアウトカムの期待値は、

$\begin{aligned} E [ y_i | x_i ] = \beta_ 0 + \beta_1 x_{1i} + \beta_2 x_{2i} \end{aligned}$

なので、説明変数 $x_1$ についての限界効果は $x_1$ で偏微分して、

$\begin{aligned} \frac{\partial}{\partial x_{1i}} E [ y_i | x_i ] = \beta_1 \end{aligned}$

となり、 線形回帰モデルでは $x_1$ の限界効果は $x_1$ 自身にも、共変量 $x_2$ にも依存せず、一定であると分かります。

ロジスティック回帰モデルの限界効果

ロジスティック回帰モデルでは、

$\begin{aligned} Pr(y_i = 1 | x_i) = \frac{1}{1+\exp \left( -\frac{\beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} }{\sigma} \right)} \end{aligned}$

なので、

$\begin{aligned} \frac{\partial}{\partial x_{1i}} Pr(y_i = 1 | x_i) = \frac{\beta_1 \exp \left( -\frac{\beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} }{\sigma} \right)}{\sigma \left( 1+\exp \left( -\frac{\beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} }{\sigma} \right) \right)^2} \end{aligned}$

となり、ロジスティック回帰モデルでは限界効果は $x_1$ にも $x_2$ にも依存していることが分かります。

ちなみにオッズ比は、

$\begin{aligned} OR = \exp \left( \frac{\beta_1}{\sigma} \right) \end{aligned}$

なので一定です。

限界効果の良いところ

オッズ比は全範囲で発生確率を適切に反映しない

下のグラフは標準ロジスティック関数です。

xが非常に小さい値の範囲では、xが1増えても発生確率はほとんど変わりませんが、中央付近だとxが1単位増えるだけで劇的に増加しています。発生確率の変化は一定でないのに、オッズ比がxの値の大小に寄らず一定というのは、本当に見たいもの（=発生確率）が反映されていない感じです。

限界効果の方がモデルに依存しにくい

オッズ比はモデルの誤差の分散に依存してしまうので、モデルに含める説明変数に依存します。これを説明変数x, zと誤差でアウトカムが決まるような仮想的なデータセットを作って考えてみましょう。

set.seed(1234)
n=10000
x <- rnorm(n,0,1)
z <- rnorm(n,0,2)

y <- NULL
for(i in 1:n){
  y[i] <- rbinom(1,1,plogis(2*x[i]+z[i]))
}

df <- data.frame(y,x,z)

このサンプルでは対数オッズに

$\begin{aligned} \log \left( \frac{p_i}{1-p_i} \right) = 2 x_i + z_i \end{aligned}$

という関係を仮定しました。これにxとzの両方を含んだモデル（full model）と、xしか含んでないモデルを当てはめます。

# 2値アウトカム
fit <- glm(y ~ x, data=df, family=binomial(link="logit"))
fit_full <- glm(y ~ x+z, data=df, family=binomial(link="logit"))

両方を含んだモデルでは、

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) 0.005211   0.028973    0.18    0.857    
x           2.016856   0.045406   44.42   <2e-16 ***
z           1.004665   0.022379   44.89   <2e-16 ***

となり、切片と各変数の係数は概ね仮定どおりに推定されています。

ここで、xしか含んでないモデルだと、

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.002755   0.022880   -0.12    0.904    
x            1.270924   0.029829   42.61   <2e-16 ***

xの係数の推定値は1.27になってしまいました（推定オッズ比にすると7.4→3.6）。

これに対して限界効果を比べると、

> library(ggeffects)
> ggpredict(fit, terms="x [all]")
# Predicted probabilities of y

    x | Predicted |       95% CI
--------------------------------
-3.40 |      0.01 | [0.01, 0.02]
-1.53 |      0.12 | [0.11, 0.14]
-0.97 |      0.23 | [0.21, 0.24]
-0.48 |      0.35 | [0.34, 0.36]
-0.01 |      0.50 | [0.49, 0.51]
 0.47 |      0.64 | [0.63, 0.65]
 0.96 |      0.77 | [0.76, 0.78]
 3.62 |      0.99 | [0.99, 0.99]

> ggpredict(fit_full, terms="x [all]")
# Predicted probabilities of y

    x | Predicted |       95% CI
--------------------------------
-3.40 |      0.00 | [0.00, 0.00]
-1.53 |      0.04 | [0.04, 0.05]
-0.97 |      0.12 | [0.11, 0.13]
-0.48 |      0.27 | [0.26, 0.29]
-0.01 |      0.49 | [0.48, 0.51]
 0.47 |      0.72 | [0.70, 0.73]
 0.96 |      0.87 | [0.86, 0.88]
 3.62 |      1.00 | [1.00, 1.00]

Adjusted for:
* z = -0.02

となり、限界効果を使う方がモデルの違いによる差があまり目立ちません。 限界効果の方が設定したモデルに影響されにくい指標だそうです。

限界効果を使うときの注意点

想定する対象者を決めなくてはならない

前述のとおり、限界効果はモデルに含まれる他の説明変数にも依存しています。そのため、「どのような患者における効果なのか」を決めなくてはなりません。 男性・女性のように変数のカテゴリー別に限界効果を求めてもいいですし、集団全体の平均値を使ってもいいです。ただし後者は、例えば「年齢が63.5歳で55%男性」がどんな人を表しているのか、よく分かりません。

ベースラインの割合と比べる必要あり

「1単位あたり2%増える」と言われても、元の割合が0.5%しかない場合と、50%の場合とでは意味合いが全く異なりますね。

ケースコントロール研究ではサンプルにおける割合が母集団の割合を反映していないので、オッズ比を使わざるをえません。

おわりに

Norton先生曰く、限界効果はリスク差みたいなもの。なるほど。
"Marginal" effectという用語に違和感を感じる人は少なくないと思います（marginalと言いながら内容はconditionalでは？）。疫学と計量経済学の分野の差でしょうか。
ggeffectsパッケージも含めようと思いましたが、想定以上に長くなってしまいました。

ランキング参加中

統計

参考資料

Norton EC, Dowd BE, Maciejewski ML. Marginal Effects-Quantifying the Effect of Changes in Risk Factors in Logistic Regression Models. JAMA. 2019 Apr 2;321(13):1304-1305.
"Marginal" という語の使い方が気になるあなたへ。

elsur.jpn.org