ロジスティック回帰

ロジスティック回帰(Logistic Regression)は二値分類および多クラス分類のための統計的モデルであり、一般化線形モデル(GLM)の枠組みにおけるベルヌーイ分布・ロジットリンクの特殊ケースとして位置づけられる。線形判別分析と異なり分布の仮定を最小限に抑えつつ、最尤推定・正則化・凸最適化の理論と直接接続する解釈可能な確率的分類モデルである。

設定

入力 $\boldsymbol{x} \in \mathbb{R}^p$、二値応答変数 $y \in \{0, 1\}$、パラメータ $\boldsymbol{\beta} \in \mathbb{R}^p$(切片を含む場合は$\boldsymbol{x}$ の第 $0$ 成分を $1$ と定義)、i.i.d. 標本 $\mathcal{D}_n = \{(\boldsymbol{x}_i, y_i)\}_{i=1}^n$ を考える。行列表記として $X \in \mathbb{R}^{n \times p}$(各行が $\boldsymbol{x}_i^\top$)を用いる。クラス $1$ の条件付き確率を

\[p_i = P(Y = 1 \mid \boldsymbol{X} = \boldsymbol{x}_i; \boldsymbol{\beta})\]

と記し、$p_i$ を $\boldsymbol{x}_i$ と $\boldsymbol{\beta}$ の関数として定式化する。

モデルの定義

シグモイド関数とロジットリンク

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

\[p_i = P(Y = 1 \mid \boldsymbol{x}_i; \boldsymbol{\beta})= \sigma(\boldsymbol{x}_i^\top \boldsymbol{\beta})= \frac{1}{1 + \exp(-\boldsymbol{x}_i^\top \boldsymbol{\beta})}\]

と定義される。$\sigma: \mathbb{R} \to (0,1)$ はシグモイド関数(ロジスティック関数)であり、以下の性質を持つ:

逆変換としてロジット関数(対数オッズ)が得られる:

\[\mathrm{logit}(p) = \log\frac{p}{1-p} = \boldsymbol{x}^\top \boldsymbol{\beta}\]

ロジット変換により $p \in (0,1)$ が $\mathbb{R}$ 全体に写され、線形予測子 $\boldsymbol{x}^\top \boldsymbol{\beta}$ との接続が実現する。オッズ(Odds)$p/(1-p)$ は「クラス $1$ になる確率のクラス $0$ に対する比」であり、$\beta_j$ は $x_j$ が $1$ 単位増加したときの対数オッズの変化量を表す。すなわち $x_j$ が $1$ 単位増加すると、オッズは $e^{\beta_j}$ 倍に変化する。

指数型分布族としての定式化

ロジスティック回帰は指数型分布族の正準形

\[p(y; \eta)= \exp\!\left(y\eta - \log(1 + e^\eta)\right)= \exp\!\left(y\eta - b(\eta)\right)\]

($\eta = \boldsymbol{x}^\top \boldsymbol{\beta}$、$b(\eta) = \log(1 + e^\eta)$、$a(\phi) = 1$)として記述できる。対数分配関数の性質より

\[\mathbb{E}[Y] = b'(\eta) = \frac{e^\eta}{1 + e^\eta} = \sigma(\eta) = p,\qquad\mathrm{Var}(Y) = b''(\eta) = \sigma(\eta)(1-\sigma(\eta)) = p(1-p)\]

が直接導かれ、正準リンク関数がロジット $g(\mu) = \log(\mu/(1-\mu))$ であることが確認できる。$b(\eta) = \log(1+e^\eta)$(Log-Sum-Exp の特殊ケース)は凸関数であり、これが後述する負の対数尤度の凸性の根拠となる。

最尤推定

尤度関数と対数尤度

$y_i \in \{0,1\}$ に対して観測の尤度は$p_i^{y_i}(1-p_i)^{1-y_i}$ であるから、対数尤度は

\[\ell(\boldsymbol{\beta})= \sum_{i=1}^n \bigl[y_i \log p_i + (1-y_i)\log(1-p_i)\bigr]= \sum_{i=1}^n \bigl[y_i \boldsymbol{x}_i^\top \boldsymbol{\beta}- \log(1 + e^{\boldsymbol{x}_i^\top \boldsymbol{\beta}})\bigr]\]

と書ける。負の対数尤度(損失関数)はバイナリ交差エントロピー損失と一致する:

\[\mathcal{L}(\boldsymbol{\beta})= -\frac{1}{n}\ell(\boldsymbol{\beta})= \frac{1}{n}\sum_{i=1}^n\bigl[-y_i\log\sigma(\boldsymbol{x}_i^\top\boldsymbol{\beta})- (1-y_i)\log(1 - \sigma(\boldsymbol{x}_i^\top\boldsymbol{\beta}))\bigr]\]

この損失関数は情報理論的には真の分布 $y_i$ とモデル予測分布 $p_i$ の間の KL ダイバージェンス(定数項を除く)に等しく、最尤推定は予測分布を真の分布に(KL 距離の意味で)最も近づける推定として解釈できる。

負の対数尤度の凸性

命題:$\mathcal{L}(\boldsymbol{\beta})$ は $\boldsymbol{\beta}$ について凸である。

証明:ヘッセ行列を計算すると

\[\nabla_{\boldsymbol{\beta}}^2 \mathcal{L}(\boldsymbol{\beta})= \frac{1}{n}\sum_{i=1}^n p_i(1-p_i)\boldsymbol{x}_i \boldsymbol{x}_i^\top= \frac{1}{n} X^\top W X\]

ここで $W = \mathrm{diag}(p_1(1-p_1), \ldots, p_n(1-p_n))$ は対角正半定値行列($p_i(1-p_i) \in (0,1/4]$ より各対角成分は正)。任意の $\boldsymbol{v} \in \mathbb{R}^p$ に対して$\boldsymbol{v}^\top \nabla^2 \mathcal{L}\boldsymbol{v} = \frac{1}{n}\|W^{1/2}X\boldsymbol{v}\|^2 \geq 0$であるから $\nabla^2\mathcal{L} \succeq 0$、すなわち $\mathcal{L}$ は凸。$\square$

$X$ が列フルランクのとき $\nabla^2\mathcal{L} \succ 0$(正定値)となり$\mathcal{L}$ は狭義凸となる——ただし後述の線形分離可能性の問題に注意を要する。凸性により局所最小解は大域最小解に一致し、勾配法・Newton 法によるアルゴリズムの大域収束が保証される。

スコア方程式と MLE の導出

対数尤度の勾配(スコア関数)は

\[\nabla_{\boldsymbol{\beta}} \ell(\boldsymbol{\beta})= \sum_{i=1}^n (y_i - p_i)\boldsymbol{x}_i= X^\top(\boldsymbol{y} - \boldsymbol{p})\]

と簡潔に表される($\boldsymbol{p} = (p_1, \ldots, p_n)^\top$)。最尤推定量 $\hat{\boldsymbol{\beta}}^{\mathrm{MLE}}$ はスコア方程式 $X^\top(\boldsymbol{y} - \boldsymbol{p}) = \boldsymbol{0}$の解であり、閉形式解は存在しないため数値的最適化が必要となる。

Newton–Raphson 法と IRLS

Newton–Raphson 法による更新則は

\[\boldsymbol{\beta}^{(t+1)}= \boldsymbol{\beta}^{(t)}- \bigl[\nabla^2 \mathcal{L}(\boldsymbol{\beta}^{(t)})\bigr]^{-1}\nabla \mathcal{L}(\boldsymbol{\beta}^{(t)})= \boldsymbol{\beta}^{(t)}+ (X^\top W^{(t)} X)^{-1} X^\top (\boldsymbol{y} - \boldsymbol{p}^{(t)})\]

と書ける。これを整理すると調整応答変数

\[\boldsymbol{z}^{(t)}= X\boldsymbol{\beta}^{(t)} + (W^{(t)})^{-1}(\boldsymbol{y} - \boldsymbol{p}^{(t)})\]

を用いた加重最小二乗問題

\[\boldsymbol{\beta}^{(t+1)}= (X^\top W^{(t)} X)^{-1} X^\top W^{(t)} \boldsymbol{z}^{(t)}= \arg\min_{\boldsymbol{\beta}}(\boldsymbol{z}^{(t)} - X\boldsymbol{\beta})^\top W^{(t)}(\boldsymbol{z}^{(t)} - X\boldsymbol{\beta})\]

と等価であり、これが GLM の一般的な解法である反復重み付き最小二乗法(IRLS)に対応する。作業重み $w_i^{(t)} = p_i^{(t)}(1-p_i^{(t)})$ は IRLS の各ステップで更新される。凸性によりこの反復は大域最適解に収束することが保証される(線形分離不可能な場合)。

線形分離可能性の問題

データが線形分離可能(Linearly Separable)である場合、すなわち $\exists \boldsymbol{\beta}^* : y_i(\boldsymbol{x}_i^\top \boldsymbol{\beta}^*) > 0$$\forall i$ を満たす $\boldsymbol{\beta}^*$ が存在する場合、MLE は存在しない。これは $\|\boldsymbol{\beta}\| \to \infty$ の方向に対数尤度が単調増加するためであり、スコア方程式の解が存在しない。幾何学的には、分離超平面の法線ベクトルの方向に$\|\boldsymbol{\beta}\| \to \infty$ とするとすべての $p_i \to y_i$ となり、対数尤度は上限 $0$ に収束するが到達しない。この問題への対処として以下が用いられる:

漸近理論

MLE の漸近正規性と有効性

線形分離可能でなく、フィッシャー情報行列

\[\mathcal{I}(\boldsymbol{\beta})= \mathbb{E}\bigl[\nabla_{\boldsymbol{\beta}}\log p(Y;\boldsymbol{x},\boldsymbol{\beta})\nabla_{\boldsymbol{\beta}}\log p(Y;\boldsymbol{x},\boldsymbol{\beta})^\top\bigr]= \mathbb{E}[p(1-p)\boldsymbol{x}\boldsymbol{x}^\top]\]

が正定値であるという正則条件のもとで、MLE は以下の漸近性質を持つ:

実用上、漸近共分散行列は$\widehat{\mathrm{Cov}}(\hat{\boldsymbol{\beta}}) = (X^\top \hat{W} X)^{-1}$($\hat{W} = \mathrm{diag}(\hat{p}_i(1-\hat{p}_i))$)で推定され、各係数の Wald 統計量

\[W_j = \frac{\hat{\beta}_j^2}{[(X^\top \hat{W} X)^{-1}]_{jj}}\xrightarrow{d} \chi^2(1)\]

が $\beta_j = 0$ の検定に用いられる。モデル全体の有意性・入れ子モデルの比較には尤度比検定(前節)が適用される。

モデル評価:逸脱度と適合度

残差逸脱度は飽和モデルとの対数尤度の差として

\[D = 2[\ell_{\mathrm{sat}} - \ell(\hat{\boldsymbol{\beta}})]= 2\sum_{i=1}^n\left[y_i\log\frac{y_i}{\hat{p}_i}+ (1-y_i)\log\frac{1-y_i}{1-\hat{p}_i}\right]\]

と定義される($y_i \in \{0,1\}$ の二値データでは$\ell_{\mathrm{sat}} = 0$)。$D$ の漸近 $\chi^2$ 近似は二値応答では一般にはよく機能せず、Hosmer–Lemeshow 検定(観測と期待頻度のグループ比較)やピアソン残差による適合度評価が補完的に用いられる。

分類モデルとしての評価指標として以下が広く用いられる。

正則化ロジスティック回帰

Ridge ロジスティック回帰($\ell_2$ 正則化)

正則化付き目的関数は

\[\mathcal{L}_\lambda(\boldsymbol{\beta})= -\frac{1}{n}\ell(\boldsymbol{\beta}) + \lambda\|\boldsymbol{\beta}\|_2^2\]

であり、$\lambda > 0$ のとき目的関数は強凸($\mu = 2\lambda$)となる。ヘッセ行列は $\nabla^2\mathcal{L}_\lambda = \frac{1}{n}X^\top WX + 2\lambda I \succ 0$であるから、線形分離可能な場合でも一意解が存在する。これはベイズ的にはガウス事前分布 $\boldsymbol{\beta} \sim \mathcal{N}(0, (2\lambda)^{-1}I)$のもとでの MAP 推定と一致する。IRLS の各ステップに Ridge ペナルティを加えた形で解かれ、更新則は

\[\boldsymbol{\beta}^{(t+1)}= (X^\top W^{(t)} X + 2n\lambda I)^{-1}(X^\top W^{(t)} X \boldsymbol{\beta}^{(t)}+ X^\top(\boldsymbol{y} - \boldsymbol{p}^{(t)}))\]

となる。漸近的には $\lambda = O(1/\sqrt{n})$ と設定することでMLE と同等の漸近有効性を保ちつつ有限標本での安定性を確保できる。

Lasso ロジスティック回帰($\ell_1$ 正則化)

$\ell_1$ 正則化は

\[\mathcal{L}_\lambda(\boldsymbol{\beta})= -\frac{1}{n}\ell(\boldsymbol{\beta}) + \lambda\|\boldsymbol{\beta}\|_1\]

と定義され、スパースな特徴選択を誘導する。$-\ell$ が凸・$\|\cdot\|_1$ が凸の和であるから目的関数は凸であり、大域最適解が存在する。$\|\cdot\|_1$ の非平滑性から IRLS は直接適用できず、近接 Newton 法(Proximal Newton Method)が用いられる:各 Newton ステップで $-\ell$ を二次近似し、近接写像(ソフト閾値処理)と組み合わせた最適化を行う。glmnet アルゴリズム(Friedman et al., 2010)は座標降下法を内側ループに用いた効率的な実装であり、$\lambda$ の正則化パス全体を $O(np\min(n,p))$ で計算する。高次元設定($p \gg n$)での変数選択一致性はアービン条件のもとで保証される。

Elastic Net ロジスティック回帰

$\ell_1$ と $\ell_2$ の混合正則化

\[\mathcal{L}_{\lambda_1,\lambda_2}(\boldsymbol{\beta})= -\frac{1}{n}\ell(\boldsymbol{\beta})+ \lambda_1\|\boldsymbol{\beta}\|_1+ \lambda_2\|\boldsymbol{\beta}\|_2^2\]

は $\ell_2$ 項による強凸性と $\ell_1$ 項によるスパース性を両立する。相関共変量が多く存在する場合に Lasso より安定した変数選択を行い、高次元バイオマーカー解析・テキスト分類等で広く用いられる。

多クラスロジスティック回帰

多項ロジスティック回帰(ソフトマックス回帰)

$K$ クラス分類($y \in \{1,\ldots,K\}$)への拡張として、クラス $k$ の係数ベクトル $\boldsymbol{\beta}_k \in \mathbb{R}^p$($k = 1,\ldots,K$)を用いて

\[P(Y = k \mid \boldsymbol{x}; B)= \frac{\exp(\boldsymbol{x}^\top \boldsymbol{\beta}_k)} {\sum_{j=1}^K \exp(\boldsymbol{x}^\top \boldsymbol{\beta}_j)},\quad k = 1, \ldots, K\]

と定義する($B = [\boldsymbol{\beta}_1, \ldots, \boldsymbol{\beta}_K] \in \mathbb{R}^{p \times K}$)。識別可能性のため $\boldsymbol{\beta}_K = \boldsymbol{0}$(参照クラス)と固定する。この場合有効パラメータは $(K-1)p$ 個であり、$K=2$ で二値ロジスティック回帰に一致する。

ソフトマックス関数$\mathrm{softmax}(\boldsymbol{\eta})_k = \exp(\eta_k)/\sum_j \exp(\eta_j)$は $\mathrm{Log\text{-}Sum\text{-}Exp}(\boldsymbol{\eta}) = \log\sum_j \exp(\eta_j)$の勾配として得られる。Log-Sum-Exp は凸関数であるから、多項ロジスティック回帰の負の対数尤度

\[\mathcal{L}(B)= -\frac{1}{n}\sum_{i=1}^n\left[\boldsymbol{x}_i^\top \boldsymbol{\beta}_{y_i}- \log\sum_{k=1}^K \exp(\boldsymbol{x}_i^\top \boldsymbol{\beta}_k)\right]\]

も $B$ について凸である(参照クラスを固定した場合は狭義凸)。スコア方程式は

\[\frac{\partial \mathcal{L}}{\partial \boldsymbol{\beta}_k}= -\frac{1}{n}\sum_{i=1}^n\bigl(\mathbf{1}[y_i = k] - \hat{p}_{ik}\bigr)\boldsymbol{x}_i= \boldsymbol{0},\quad k = 1, \ldots, K-1\]

と書かれ($\hat{p}_{ik} = P(Y=k\mid\boldsymbol{x}_i;\hat{B})$)、IRLS または L-BFGS 等の準 Newton 法で数値的に解かれる。

One-vs-Rest と One-vs-One

多クラス分類への別アプローチとして、二値ロジスティック回帰を繰り返し適用する戦略がある。

確率を一貫して推定するには多項ロジスティック回帰(ソフトマックス回帰)が理論的に適切であり、OvR・OvO は近似的手法として位置づけられる。

ロジスティック回帰と線形判別分析の比較

二値分類においてロジスティック回帰と線形判別分析(LDA)は、どちらも線形決定境界 $\boldsymbol{x}^\top \boldsymbol{\beta} = 0$ を生成するが、推定のアプローチが異なる。

観点ロジスティック回帰線形判別分析(LDA)
モデル化の対象$P(Y \mid \boldsymbol{X})$(判別的モデル)$P(\boldsymbol{X} \mid Y)$ と $P(Y)$(生成モデル)
分布の仮定$P(\boldsymbol{X})$ に仮定なし$\boldsymbol{X} \mid Y=k \sim \mathcal{N}(\boldsymbol{\mu}_k, \Sigma)$
推定法条件付き尤度の最大化同時尤度の最大化(閉形式)
小標本・仮定満足時LDA より劣る傾向漸近有効(正規性仮定が真のとき)
大標本・仮定違反時ロバストで優れた性能正規性・等共分散仮定違反で劣化
パラメータ解釈対数オッズ比として直接解釈可能分布パラメータ経由の間接解釈
高次元・正則化$\ell_1$/$\ell_2$ 正則化が自然に適用可能正則化 LDA(RLDA)が必要

正規性仮定が成立するとき LDA はロジスティック回帰より漸近的に有効(少ない標本で同精度)であるが、仮定違反に対してロジスティック回帰はよりロバストである(Efron, 1975 による漸近相対効率の分析)。実用上は分布の仮定が疑わしい場合にロジスティック回帰が選好される。

ロジスティック回帰の確率的解釈と潜在変数モデル

ロジスティック回帰は潜在変数モデルとして解釈できる。潜在変数 $Z = \boldsymbol{x}^\top \boldsymbol{\beta} + \varepsilon$($\varepsilon$ は標準ロジスティック分布 $F(\varepsilon) = \sigma(\varepsilon)$ に従う)として、

\[Y = \mathbf{1}[Z > 0]\impliesP(Y=1 \mid \boldsymbol{x})= P(Z > 0)= P(\varepsilon > -\boldsymbol{x}^\top\boldsymbol{\beta})= \sigma(\boldsymbol{x}^\top\boldsymbol{\beta})\]

が成立する。$\varepsilon$ を標準正規分布 $\mathcal{N}(0,1)$ に置き換えるとプロビットモデルが得られる:$P(Y=1\mid\boldsymbol{x}) = \Phi(\boldsymbol{x}^\top\boldsymbol{\beta})$。ロジスティック分布と正規分布は尾部の厚さが異なり(ロジスティック分布の方が裾が重い)、両モデルの予測はほぼ一致するが係数のスケールが異なる(係数比の近似 $\beta_{\mathrm{logistic}} \approx (\pi/\sqrt{3})\, \beta_{\mathrm{probit}}$)。

ベイズロジスティック回帰においては、Polya-Gamma 補助変数(Polson–Scott–Windle, 2013)を用いたデータ拡張により、事後分布からの効率的なギブスサンプリングが実現する。具体的には $\omega_i \sim \mathrm{PG}(1, \boldsymbol{x}_i^\top\boldsymbol{\beta})$ を補助変数として導入することで、条件付き事後分布がガウスとなり、共役更新が可能となる。

まとめ

ロジスティック回帰はシグモイド関数によって線形予測子を条件付き確率に変換する判別的確率モデルであり、ベルヌーイ・ロジットリンクの GLM として指数型分布族の正準形に位置づけられる。負の対数尤度(バイナリ交差エントロピー)の凸性が大域的な最適化を保証し、IRLS(Newton–Raphson 法)による効率的な計算を可能にする。MLE は漸近正規性・有効性・一致性を持ち、Wald 検定・尤度比検定による統計的推測の枠組みを提供する。線形分離可能性の問題は $\ell_2$ 正則化により解決され、$\ell_1$ 正則化は高次元設定でのスパース変数選択を実現する。多項ロジスティック回帰はソフトマックス関数を通じて$K$ クラス問題へと自然に拡張され、ニューラルネットワークの最終層の理論的基礎をなす。LDA との比較では分布の仮定に関するロバスト性と漸近有効性のトレードオフが存在し、潜在変数モデルとしての解釈はプロビットモデル・ベイズ推定との接続を与える。

Mathematics is the language with which God has written the universe.





















指数型分布族 主成分分析(PCA) 因子分析 判別分析(LDA/QDA) 正準相関分析(CCA)