多項ロジスティック回帰

多項ロジスティック回帰の基礎となる数学モデルを記述する方法は複数あり、それぞれ同等の表現が存在します。そのため、異なる文献におけるこのテーマの扱い方を比較することは困難です。ロジスティック回帰に関する本稿では、単純ロジスティック回帰の同等の定式化をいくつか紹介しており、これらの多くは多項ロジットモデルにも類似の表現が存在します。

他の多くの統計的分類手法と同様に、これらの手法の根底にある考え方は、与えられた観測値の説明変数（特徴）と内積を用いて線形結合された一連の重みからスコアを構築する線形予測関数を構築することである。

${\displaystyle \operatorname {score} (\mathbf {X} _{i},k)={\boldsymbol {\beta }}_{k}\cdot \mathbf {X} _{i},}$

ここで、X _iは観測iを説明する説明変数のベクトル、β _kは結果kに対応する重み（または回帰係数）のベクトル、score( X _i , k ) は観測i をカテゴリkに割り当てることに関連付けられたスコアです。離散選択理論では、観測は人を表し、結果は選択を表すため、スコアは人i が結果kを選択することに関連付けられた効用とみなされます。予測された結果は、最も高いスコアを持つものです。

多項ロジットモデルと、同じ基本設定を持つ他の多くの手法、モデル、アルゴリズム（パーセプトロンアルゴリズム、サポートベクターマシン、線形判別分析など）との違いは、最適な重み／係数を決定（学習）する手順と、スコアの解釈方法にあります。特に、多項ロジットモデルでは、スコアを直接確率値に変換できます。これは、観測の測定された特性に基づいて、観測iが結果kを選択する確率を示します。これにより、特定の多項ロジットモデルの予測を、それぞれに誤差が生じる可能性のある複数の予測を含む、より大きな手順に組み込むための原理に基づいた方法が提供されます。このような予測を組み合わせる手段がない場合、誤差は増幅する傾向があります。例えば、大規模な予測モデルを、一連のサブモデルに分割し、あるサブモデルの予測結果を別のサブモデルの入力として使用し、その予測結果をさらに別のサブモデルの入力として使用する、といったように分割していく場合を考えてみましょう。各サブモデルの予測精度が90%で、5つのサブモデルが直列に接続されている場合、全体のモデルの精度は0.95 = 59%にしかなりません^。各サブモデルの精度が80%の場合、全体の精度は0.85 = 33%に低下します^。この問題は誤差伝播として知られており、通常多数の部分から構成される現実世界の予測モデルでは深刻な問題となります。この問題を軽減する一つの方法として、単一の最適な予測を行うのではなく、起こりうる各結果の確率を予測することが挙げられます。

基本的な設定はロジスティック回帰と同じですが、唯一の違いは従属変数が二値変数ではなくカテゴリ変数である点です。つまり、結果は2つではなくK個存在することになります。以下の説明はやや簡略化されています。詳細については、ロジスティック回帰に関する記事を参照してください。

多項ロジットモデルを導出するには、K個の可能な結果に対して、K個の独立した二項ロジスティック回帰モデルを実行することを想定できます。このモデルでは、1つの結果を「ピボット」として選択し、残りのK  - 1個の結果をピボット結果に対して個別に回帰分析します。結果K（最後の結果）をピボットとして選択した 場合、 K - 1個の回帰方程式は次のようになります。  

${\displaystyle \ln {\frac {\Pr(Y_{i}=k)}{\Pr(Y_{i}=K)}}\,=\,{\boldsymbol {\beta }}_{k}\cdot \mathbf {X} _{i},\;\;\;\;\;\;1\leq k<K}$。

この定式化は、構成データ分析で一般的に使用される加法対数比変換としても知られています。他のアプリケーションでは、「相対リスク」と呼ばれています。^{[ 7 ]}

両辺を指数化して確率を求めると、次の式が得られます。

${\displaystyle \Pr(Y_{i}=k)\,=\,{\Pr(Y_{i}=K)}\;e^{{\boldsymbol {\beta }}_{k}\cdot \mathbf {X} _{i}},\;\;\;\;\;\;1\leq k<K}$

K個の確率の合計が1になるという事実を用いると、次のことがわかる。

${\displaystyle {\begin{aligned}\Pr(Y_{i}=K)={}&1-\sum _{j=1}^{K-1}\Pr(Y_{i}=j)\\={}&1-\sum _{j=1}^{K-1}{\Pr(Y_{i}=K)}\;e^{{\boldsymbol {\beta }}_{j}\cdot \mathbf {X} _{i}}\;\;\Rightarrow \;\;\Pr(Y_{i}=K)\\={}&{\frac {1}{1+\sum _{j=1}^{K-1}e^{{\boldsymbol {\beta }}_{j}\cdot \mathbf {X} _{i}}}}.\end{aligned}}}$

これを利用して、他の確率を求めることができます。

${\displaystyle \Pr(Y_{i}=k)={\frac {e^{{\boldsymbol {\beta }}_{k}\cdot \mathbf {X} _{i}}}{1+\sum _{j=1}^{K-1}e^{{\boldsymbol {\beta }}_{j}\cdot \mathbf {X} _{i}}}},\;\;\;\;\;\;1\leq k<K}$。

複数の回帰分析を実行するという事実は、なぜこのモデルが上述の無関係な選択肢の独立性という仮定に基づいているのかを明らかにしている。

The unknown parameters in each vector β_k are typically jointly estimated by maximum a posteriori (MAP) estimation, which is an extension of maximum likelihood using regularization of the weights to prevent pathological solutions (usually a squared regularizing function, which is equivalent to placing a zero-mean Gaussianprior distribution on the weights, but other distributions are also possible). The solution is typically found using an iterative procedure such as generalized iterative scaling,^[8]iteratively reweighted least squares (IRLS),^[9] by means of gradient-based optimization algorithms such as L-BFGS,^[4] or by specialized coordinate descent algorithms.^[10]

The formulation of binary logistic regression as a log-linear model can be directly extended to multi-way regression. That is, we model the logarithm of the probability of seeing a given output using the linear predictor as well as an additional normalization factor, the logarithm of the partition function:

${\displaystyle \ln \Pr(Y_{i}=k)={\boldsymbol {\beta }}_{k}\cdot \mathbf {X} _{i}-\ln Z,\;\;\;\;\;\;1\leq k\leq K.}$

As in the binary case, we need an extra term ${\displaystyle -\ln Z}$ to ensure that the whole set of probabilities forms a probability distribution, i.e. so that they all sum to one:

${\displaystyle \sum _{k=1}^{K}\Pr(Y_{i}=k)=1}$

The reason why we need to add a term to ensure normalization, rather than multiply as is usual, is because we have taken the logarithm of the probabilities. Exponentiating both sides turns the additive term into a multiplicative factor, so that the probability is just the Gibbs measure:

${\displaystyle \Pr(Y_{i}=k)={\frac {1}{Z}}e^{{\boldsymbol {\beta }}_{k}\cdot \mathbf {X} _{i}},\;\;\;\;\;\;1\leq k\leq K.}$

The quantity Z is called the partition function for the distribution. We can compute the value of the partition function by applying the above constraint that requires all probabilities to sum to 1:

${\displaystyle 1=\sum _{k=1}^{K}\Pr(Y_{i}=k)\;=\;\sum _{k=1}^{K}{\frac {1}{Z}}e^{{\boldsymbol {\beta }}_{k}\cdot \mathbf {X} _{i}}\;=\;{\frac {1}{Z}}\sum _{k=1}^{K}e^{{\boldsymbol {\beta }}_{k}\cdot \mathbf {X} _{i}}.}$

Therefore

${\displaystyle Z=\sum _{k=1}^{K}e^{{\boldsymbol {\beta }}_{k}\cdot \mathbf {X} _{i}}.}$

Note that this factor is "constant" in the sense that it is not a function of Y_i, which is the variable over which the probability distribution is defined. However, it is definitely not constant with respect to the explanatory variables, or crucially, with respect to the unknown regression coefficients β_k, which we will need to determine through some sort of optimization procedure.

The resulting equations for the probabilities are

${\displaystyle \Pr(Y_{i}=k)={\frac {e^{{\boldsymbol {\beta }}_{k}\cdot \mathbf {X} _{i}}}{\sum _{j=1}^{K}e^{{\boldsymbol {\beta }}_{j}\cdot \mathbf {X} _{i}}}},\;\;\;\;\;\;1\leq k\leq K.}$

The following function:

${\displaystyle \operatorname {softmax} (k,s_{1},\ldots ,s_{K})={\frac {e^{s_{k}}}{\sum _{j=1}^{K}e^{s_{j}}}}}$

これはソフトマックス関数と呼ばれます。その理由は、値を指数化する効果が${\displaystyle s_{1},\ldots ,s_{K}}$両者の違いを誇張することである。その結果、${\displaystyle \operatorname {softmax} (k,s_{1},\ldots ,s_{K})}$0に近い値を返す${\displaystyle s_{k}}$はすべての値の最大値よりもかなり小さく、最大値に適用すると、次に大きい値に極めて近い場合を除き、1に近い値を返します。したがって、ソフトマックス関数は、滑らかな関数（簡単に微分できるなど）として振る舞い、指示関数を近似する加重平均を構築するために使用できます。

${\displaystyle f(k)={\begin{cases}1&{\textrm {if}}\;k=\operatorname {\arg \max } _{j}s_{j},\\0&{\textrm {otherwise}}.\end{cases}}}$

したがって、確率方程式は次のように書くことができます。

${\displaystyle \Pr(Y_{i}=k)=\operatorname {softmax} (k,{\boldsymbol {\beta }}_{1}\cdot \mathbf {X} _{i},\ldots ,{\boldsymbol {\beta }}_{K}\cdot \mathbf {X} _{i})}$

したがって、ソフトマックス関数は、二項ロジスティック回帰におけるロジスティック関数と同等の役割を果たす。

すべての${\displaystyle {\boldsymbol {\beta }}_{k}}$係数ベクトルは一意に識別可能です。これは、すべての確率の合計が 1 になる必要があるため、残りのすべてがわかれば、そのうちの 1 つが完全に決定されるためです。結果として、係数ベクトルは 1 つしか存在しません。${\displaystyle K-1}$個別に指定可能な確率、したがって${\displaystyle K-1}$係数のベクトルはそれぞれ識別可能である。これを理解する一つの方法は、すべての係数ベクトルに定数ベクトルを加えると、方程式が同一になることに注目することである。

${\displaystyle {\begin{aligned}{\frac {e^{({\boldsymbol {\beta }}_{k}+\mathbf {C} )\cdot \mathbf {X} _{i}}}{\sum _{j=1}^{K}e^{({\boldsymbol {\beta }}_{j}+\mathbf {C} )\cdot \mathbf {X} _{i}}}}&={\frac {e^{{\boldsymbol {\beta }}_{k}\cdot \mathbf {X} _{i}}e^{\mathbf {C} \cdot \mathbf {X} _{i}}}{\sum _{j=1}^{K}e^{{\boldsymbol {\beta }}_{j}\cdot \mathbf {X} _{i}}e^{\mathbf {C} \cdot \mathbf {X} _{i}}}}\\&={\frac {e^{\mathbf {C} \cdot \mathbf {X} _{i}}e^{{\boldsymbol {\beta }}_{k}\cdot \mathbf {X} _{i}}}{e^{\mathbf {C} \cdot \mathbf {X} _{i}}\sum _{j=1}^{K}e^{{\boldsymbol {\beta }}_{j}\cdot \mathbf {X} _{i}}}}\\&={\frac {e^{{\boldsymbol {\beta }}_{k}\cdot \mathbf {X} _{i}}}{\sum _{j=1}^{K}e^{{\boldsymbol {\beta }}_{j}\cdot \mathbf {X} _{i}}}}\end{aligned}}}$

その結果、設定するのが一般的である${\displaystyle \mathbf {C} =-{\boldsymbol {\beta }}_{K}}$（または、他の係数ベクトルのいずれか）。基本的に、定数を設定して、ベクトルの1つが${\displaystyle {\boldsymbol {0}}}$そして、他のすべてのベクトルは、それらのベクトルと選択したベクトルとの差に変換されます。これは、K 個の選択肢のうちの 1 つを中心に「ピボット」し、他のK  −  1 個の選択肢が、中心となる選択肢と比較してどれだけ優れているか、あるいは劣っているかを調べることに相当します。数学的には、係数は次のように変換されます。

${\displaystyle {\begin{aligned}{\boldsymbol {\beta }}'_{k}&={\boldsymbol {\beta }}_{k}-{\boldsymbol {\beta }}_{K},\;\;\;\;1\leq k<K,\\{\boldsymbol {\beta }}'_{K}&=0.\end{aligned}}}$

これにより、以下の式が得られます。

${\displaystyle \Pr(Y_{i}=k)={\frac {e^{{\boldsymbol {\beta }}'_{k}\cdot \mathbf {X} _{i}}}{1+\sum _{j=1}^{K-1}e^{{\boldsymbol {\beta }}'_{j}\cdot \mathbf {X} _{i}}}},\;\;\;\;\;\;1\leq k\leq K}$

回帰係数にプライム記号が付いている点を除けば、これはK  -  1個の独立した2方向回帰という観点から、上述のモデルの形式と全く同じである。

多項ロジスティック回帰は、二項ロジスティック回帰で説明した二方向潜在変数モデルに倣って、潜在変数モデルとして定式化することも可能です。この定式化は離散選択モデルの理論において一般的であり、多項ロジスティック回帰を関連する多項プロビットモデルと比較しやすくするだけでなく、より複雑なモデルへの拡張も容易にします。

各データポイントiと可能な結果k  =  1,2,..., Kに対して、次のような分布に従う連続潜在変数Y _i,k ^* (つまり、観測されない確率変数) が存在すると想像してください。

${\displaystyle Y_{i,k}^{\ast }={\boldsymbol {\beta }}_{k}\cdot \mathbf {X} _{i}+\varepsilon _{k}\;\;\;\;,\;\;k\leq K}$

どこ${\displaystyle \varepsilon _{k}\sim \operatorname {EV} _{1}(0,1),}$つまり、標準的なタイプ1の極値分布です。

この潜在変数は、データポイントiが結果kを選択した場合の効用と考えることができます。ただし、実際に得られる効用の量にはある程度のランダム性があり、これは選択に影響を与える他のモデル化されていない要因を考慮に入れています。実際の変数の値は${\displaystyle Y_{i}}$は、これらの潜在変数から非ランダムな方法で決定されます（つまり、ランダム性は観測された結果から潜在変数に移されています）。結果kは、関連する効用（値）が${\displaystyle Y_{i,k}^{\ast }}$）は他のすべての選択肢の効用よりも大きい、つまり結果kに関連付けられた効用がすべての効用の最大値である場合です。潜在変数は連続であるため、2 つの値がまったく同じになる確率は 0 なので、このシナリオは無視します。つまり、次のようになります。

${\displaystyle {\begin{aligned}\Pr(Y_{i}=1)&=\Pr(Y_{i,1}^{\ast }>Y_{i,2}^{\ast }{\text{ and }}Y_{i,1}^{\ast }>Y_{i,3}^{\ast }{\text{ and }}\cdots {\text{ and }}Y_{i,1}^{\ast }>Y_{i,K}^{\ast })\\\Pr(Y_{i}=2)&=\Pr(Y_{i,2}^{\ast }>Y_{i,1}^{\ast }{\text{ and }}Y_{i,2}^{\ast }>Y_{i,3}^{\ast }{\text{ and }}\cdots {\text{ and }}Y_{i,2}^{\ast }>Y_{i,K}^{\ast })\\&\,\,\,\vdots \\\Pr(Y_{i}=K)&=\Pr(Y_{i,K}^{\ast }>Y_{i,1}^{\ast }{\text{ and }}Y_{i,K}^{\ast }>Y_{i,2}^{\ast }{\text{ and }}\cdots {\text{ and }}Y_{i,K}^{\ast }>Y_{i,K-1}^{\ast })\\\end{aligned}}}$

あるいは同等に：

${\displaystyle \Pr(Y_{i}=k)\;=\;\Pr(\max(Y_{i,1}^{\ast },Y_{i,2}^{\ast },\ldots ,Y_{i,K}^{\ast })=Y_{i,k}^{\ast })\;\;\;\;,\;\;k\leq K}$

最初の式を詳しく見てみましょう。これは次のように書くことができます。

${\displaystyle {\begin{aligned}\Pr(Y_{i}=1)&=\Pr(Y_{i,1}^{\ast }>Y_{i,k}^{\ast }\ \forall \ k=2,\ldots ,K)\\&=\Pr(Y_{i,1}^{\ast }-Y_{i,k}^{\ast }>0\ \forall \ k=2,\ldots ,K)\\&=\Pr({\boldsymbol {\beta }}_{1}\cdot \mathbf {X} _{i}+\varepsilon _{1}-({\boldsymbol {\beta }}_{k}\cdot \mathbf {X} _{i}+\varepsilon _{k})>0\ \forall \ k=2,\ldots ,K)\\&=\Pr(({\boldsymbol {\beta }}_{1}-{\boldsymbol {\beta }}_{k})\cdot \mathbf {X} _{i}>\varepsilon _{k}-\varepsilon _{1}\ \forall \ k=2,\ldots ,K)\end{aligned}}}$

ここで留意すべき点がいくつかあります。

1. 一般的に、もし${\displaystyle X\sim \operatorname {EV} _{1}(a,b)}$そして${\displaystyle Y\sim \operatorname {EV} _{1}(a,b)}$それから${\displaystyle X-Y\sim \operatorname {Logistic} (0,b).}$つまり、独立で同一分布に従う2つの極値分布変数の差はロジスティック分布に従い、最初のパラメータは重要ではありません。これは、最初のパラメータが位置パラメータ、つまり平均値を一定量だけシフトさせるものであり、2つの値が両方とも同じ量だけシフトされた場合、その差は同じままであることから理解できます。これは、特定の選択の確率の根底にあるすべての関係式がロジスティック分布に関係していることを意味し、一見恣意的に思えた極値分布という最初の選択が、いくらか理解しやすくなります。
2. 極値分布またはロジスティック分布の 2 番目のパラメータはスケール パラメータであり、${\displaystyle X\sim \operatorname {Logistic} (0,1)}$それから${\displaystyle bX\sim \operatorname {Logistic} (0,b).}$これは、スケール1の代わりに任意のスケールパラメータを持つ誤差変数を使用した場合の影響は、すべての回帰ベクトルに同じスケールを乗じるだけで簡単に補償できることを意味します。前の点と合わせて考えると、誤差変数に標準的な極値分布（位置0、スケール1）を使用しても、任意の極値分布を使用した場合と比べて一般性が損なわれることはないことがわかります。実際、より一般的な分布を使用した場合、モデルは識別不可能（最適な係数の単一のセットが存在しない）になります。
3. 回帰係数ベクトルの差分のみが使用されるため、すべての係数ベクトルに任意の定数を加えてもモデルには影響がありません。つまり、対数線形モデルと同様に、係数ベクトルのうちK  -  1個のみが識別可能であり、最後の1つは任意の値（例えば0）に設定できます。

実際に上記の確率の値を求めるのはやや難しく、一連の値の中から特定の順序統計量（最初の値、つまり最大値）を計算する問題となります。しかし、得られる式は上記の定式化と同じであり、両者は等価であることが示せます。