はじめに

この記事では、最尤推定法（MLE）を用いたロジット係数の推定と、これらの係数をオッズ比に変換する方法について詳しく説明し、ロジットモデルの評価と検証についても紹介します。

また、モデルの推定・解釈・評価をRを使用して実践的なデモンストレーションを行います。

ロジット係数の推定と解釈

この章では、最尤推定法（MLE）を用いたロジット係数の推定と、これらの係数をオッズ比に変換する方法について説明します。

最尤推定法

ロジットモデルでは、2値の結果変数 Y と一連の予測変数 X_1, X_2, \dots, X_p の関係は、オッズ比率の自然対数であるロジット関数で表されます。

係数 \beta_0, \beta_1, \dots, \beta_p を推定するには、最尤法（MLE）を用います。ロジットモデルの尤度関数は以下のようになります。

最尤推定量は、尤度関数を最大にする係数の値となります。これらの推定量を求めるために、通常、Newton-Raphson法や反復重み付き最小二乗法（IRLS）などの反復的な数値最適化アルゴリズムが用いられます。

オッズ比と解釈

ロジット係数を解釈するために、しばしばオッズ比に変換します。オッズ比は、予測変数の2つの異なる値に対する応答変数が1であるオッズの比率です。予測変数 X_j の1単位増加に対するオッズ比は以下のようになります。

1より大きいオッズ比率は、予測変数の1単位増加により結果がより起こりやすくなることを示し、1より小さいオッズ比率は、結果がより起こりにくくなることを示します。オッズ比率が1の場合、予測変数が結果変数に影響を与えないことを示します。

オッズ比率の解釈を理解するために、次の例を考えてみます。年齢と体重指数（BMI）に基づいて、糖尿病を持つ可能性を推定するロジットモデルがあるとします。推定されたロジットモデルの係数が \beta_1 = 0.05（年齢）および \beta_2 = 0.15（BMI）であるとします。

年齢のオッズ比率は e^{0.05} \approx 1.05 であり、1年の年齢が増えるごとに、糖尿病の発症率が約5％増加することを示します。BMIのオッズ比率は e^{0.15} \approx 1.16 であり、BMIが1単位増加するごとに、糖尿病の発症率が約16％増加することを示します。

モデル評価と検証

ロジットモデルを推定した後、パフォーマンスを評価し、妥当性を評価することが重要です。この章では、適合度の良し悪しを示す尺度について説明し、モデルの仮定や限界を検討します。

適合度の良し悪しを示す尺度

ロジットモデルの適合度を評価するためには、尤度比検定、アカイケ情報量基準（AIC）、ベイズ情報量基準（BIC）などの尺度が使用されます。これらの尺度は、異なるモデルの適合度を比較し、予測変数を追加または削除することでモデルを改善するかどうかを判断するのに役立ちます。

尤度比検定

尤度比検定は、1つのモデルが他のモデルの部分集合である、2つの入れ子のモデルの適合度を比較します。検定統計量は次のように与えられます。

ここで、L_0 および L_1 は、それぞれヌルモデルと代替モデルの尤度です。検定統計量は、2つのモデル間のパラメータ数の差に等しい自由度を持つカイ二乗分布に従います。

アカイケ情報量基準（AIC）

AICは、適合度とモデルの複雑さをバランスするモデル適合度の尺度です。AICの値が小さいほど、適合度が良いモデルです。AICは次のように与えられます。

ここで、L はモデルの尤度、k は推定されたパラメータ数です。

ベイズ情報量基準（BIC）

AICと同様に、BICも適合度とモデルの複雑さをバランスする尺度ですが、パラメータの追加に対してより強いペナルティを課します。BICの値が小さいほど、適合度が良いモデルです。BICは次のように与えられます。

ここで、n はサンプルサイズです。

疑似R2

疑似 R^2 値は、McFaddenの R^2 など、線形回帰の R^2 と比較できるモデル適合度の代替尺度を提供します。McFaddenの R^2 は次のように与えられます。

ここで、L_0 はヌルモデル（切片のみ）の尤度であり、L_1 は推定されたモデルの尤度です。

モデルの仮定と限界

ロジットモデルには、解釈結果に考慮する必要がある仮定や限界があります。

• ロジットの線形性
  ロジットモデルは、目的変数の確率のロジット関数が予測変数と線形関係にあることを仮定しています。この仮定は全てのケースで成立しない場合があり、予測変数を変換したり、相互作用項を含めたりする必要がある場合があります。

• 観測値の独立性
  ロジットモデルは、観測値が独立であると仮定しています。ロングデータやクラスターデータなど、観測値に相関がある場合は、混合効果モデルや一般化推定方程式（GEE）などの特殊な手法を考慮する必要があります。

• 完全な分離のないこと
  ロジットモデルは、予測変数の線形結合による目的変数の完全分離がないことを仮定しています。完全分離があると、ロジット係数の推定値が無限になる可能性があります。

• 大規模なサンプルサイズ
  ロジットモデルは、最大尤度推定値の妥当性と標準誤差の推定に大規模なサンプルサイズが必要です。サンプルサイズが小さい場合、推定値に偏りが生じ、信頼区間が不正確になる可能性があります。そのような場合には、ペナルティ付き最尤推定法やベイズ法などの代替的な推定法がより適している場合があります。

• 多重共線性
  ロジットモデルは、他の回帰モデルと同様に、予測変数間の多重共線性に敏感に反応することがあります。多重共線性があると、不安定な推定値、膨張した標準誤差、係数の解釈が困難になる可能性があります。多重共線性をチェックし、高度に相関する予測変数を削除または結合するか、主成分分析（PCA）のような次元削減技術を使用することが重要です。

R によるロジットモデルの推定と解釈

ロジット係数の推定と解釈、そしてRを使ったモデルの評価の例を示します。

データの準備

まず、必要なライブラリを読み込んで、シミュレートされたデータセットを作成します。

library(dplyr)
library(ggplot2)
library(caret)

set.seed(123)

n <- 1000
age <- rnorm(n, mean = 45, sd = 10)
bmi <- rnorm(n, mean = 25, sd = 5)
probability <- exp(0.05 * age + 0.15 * bmi) / (1 + exp(0.05 * age + 0.15 * bmi))
has_diabetes <- rbinom(n, size = 1, prob = probability)

data <- data.frame(has_diabetes, age, bmi)

ロジットモデルの推定

次に、glm()関数を使用してロジットモデルを推定します。

logit_model <- glm(has_diabetes ~ age + bmi, data = data, family = binomial(link = "logit"))
summary(logit_model)

Call:
glm(formula = has_diabetes ~ age + bmi, family = binomial(link = "logit"),
    data = data)

Deviance Residuals:
    Min       1Q   Median       3Q      Max
-3.1832   0.0421   0.0646   0.1021   0.6158

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.13537    2.52900  -0.844   0.3985
age          0.11552    0.04915   2.350   0.0188 *
bmi          0.12255    0.09559   1.282   0.1998
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 62.958  on 999  degrees of freedom
Residual deviance: 54.338  on 997  degrees of freedom
AIC: 60.338

Number of Fisher Scoring iterations: 9

ロジット係数の解釈

オッズ比を計算することで、ロジット係数を解釈することができます。

exp(coef(logit_model))

(Intercept)         age         bmi
  0.1182008   1.1224517   1.1303729

モデルの評価

モデルの適合度を評価するために、様々な適合度指標を使用することができます。

尤度比検定

null_model <- glm(has_diabetes ~ 1, data = data, family = binomial(link = "logit"))
anova(null_model, logit_model, test = "Chisq")

Analysis of Deviance Table

Model 1: has_diabetes ~ 1
Model 2: has_diabetes ~ age + bmi
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1       999     62.958
2       997     54.338  2     8.62  0.01343 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

AIC と BIC

AIC(logit_model)
BIC(logit_model)

[1] 60.33809
[1] 75.06135

McFadden の擬似R2

1 - logLik(logit_model) / logLik(null_model)

'log Lik.' 0.1369171 (df=3)