累積ロジットモデル（cumulative logit model）

比例オッズモデル（proportional odds model） や 累積ロジットモデル（cumulative logit model） と呼ばれる

順序ロジットモデル（ordered logit model） とも呼ばれる様子（Ordered logit - Wikipedia）

潜在連続変数モデル（Latent Variable Formulation）¶

順序ロジットでは、 $y \in \{1,2,\dots,K\}$ という順序カテゴリは、背後に潜在連続変数

y^* = \mathbf{x}^\top \beta + \varepsilon ,\quad \varepsilon \sim \text{Logistic}(0,1)

(1)

が存在し、それを閾値（cutpoint） $\tau_k$ で区切った結果であると考える：

y = k \iff \tau_{k-1} < y^* \le \tau_k

(2)

比例オッズモデル（proportional odds model）¶

この仮定のもとで、累積確率は次のように表される：

P(y \le k) = \frac{1} {1+\exp\left[-(\tau_k - \mathbf{x}^\top \beta)\right]} \quad (k=1,\dots,K-1)

(5)

各カテゴリに属する確率 $P(y = k)$ は、累積確率の差で表現する。

P(y = k) = P(y \le k) -P(y \le k-1)

(6)

つまり

P(y = k) = \frac{1} {1+\exp\left[-(\tau_k - \mathbf{x}^\top \beta)\right]} - \frac{1} {1+\exp\left[-(\tau_{k-1} - \mathbf{x}^\top \beta)\right]}

(7)

と表される。

端のカテゴリの扱い¶

両端の閾値は $-\infty, \infty$ と仮定する。

-\infty \equiv \tau_0 \leq \tau_1 \leq \ldots \leq \tau_{K-1} \leq \tau_K \equiv \infty

(8)

よって、両端の確率は $0, 1$ と仮定する

P(y \le 0) = 0,\quad P(y \le K) = 1

(9)

こうすることで、

最小カテゴリ： $P(y=1)=P(y \le 1)$
最大カテゴリ： $P(y=K)=1 - P(y \le K-1)$

と書ける。

なお、累積確率の両辺のロジットをとると

\log \frac{P(y \le k)} {P(y > k)} = \tau_k - \mathbf{x}^\top \beta \quad (k=1,\dots,K-1)

(10)

と表すこともできる

比例オッズ仮定（Proportional Odds Assumption）¶

比例オッズモデル（proportional odds model）と累積ロジットモデル（cumulative logit model）を区別しない文献もあるが、区別する文献においては、比例オッズモデルは

切片（ $\tau_k$ ）は $k$ ごとに異なる
傾き（ $\beta$ ）はすべて共通

という制約（parallel slopes constraint）を課すモデルであるとされる。

そして累積ロジットモデルではカテゴリごとに異なる傾き $\beta_k$ を持つ

\log \frac{P(y \le k)} {P(y > k)} = \mathbf{x}^\top \beta_k \quad (k=1,\dots,K-1)

(12)

実装¶

サンプルデータ¶

gpa：GPA。 $[0,4]$ のfloat
pared：少なくとも片方の親は大学院に進学している。二値変数。
public：学生が現在在籍している大学がpublicかprivateか
apply： unlikely < somewhat likely < very likely の順序変数

import pandas as pd
url = "https://stats.idre.ucla.edu/stat/data/ologit.dta"
df = pd.read_stata(url)

y = df['apply']
X = df[['pared', 'public', 'gpa']]

df.head(5)

Loading...

statsmodelsによる実装（最尤推定）¶

import numpy as np
import pandas as pd
import scipy.stats as stats

from statsmodels.miscmodels.ordinal_model import OrderedModel

model = OrderedModel(y, X, distr='logit')
result = model.fit(method='bfgs')
result.summary()

Optimization terminated successfully.
         Current function value: 0.896281
         Iterations: 22
         Function evaluations: 24
         Gradient evaluations: 24

Loading...

PyMCによる実装（MCMC）¶

PyMCのOrderedLogisticでは予測子（predictor）をeta $\eta = x^\top \beta$ 、閾値 $\tau$ をcutpointsで表す。

参考：

import numpy as np
import pymc as pm
import pytensor.tensor as pt
from statsmodels.miscmodels.ordinal_model import OrderedModel

N, M = X.shape
y_num = y.cat.codes

with pm.Model() as model:
    beta = pm.Normal("beta", mu=0, sigma=1, shape=M)
    eta = X.to_numpy() @ beta
    cutpoints = pm.Normal("cutpoints", mu=[-1,1], sigma=10, shape=2,
                          transform=pm.distributions.transforms.ordered)
    y_ = pm.OrderedLogistic("y", cutpoints=cutpoints, eta=eta, observed=y_num)
    idata = pm.sample()

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [beta, cutpoints]

Loading...

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 5 seconds.

# EAP（事後平均）推定量
beta_EAP = idata.posterior["beta"].mean(dim=("chain", "draw"))
pd.DataFrame({
    "exog": X.columns,
    "EAP": beta_EAP.to_numpy().round(4),
})

Loading...

潜在連続変数モデル（Latent Variable Formulation）¶

比例オッズモデル（proportional odds model）¶

端のカテゴリの扱い¶

比例オッズ仮定（Proportional Odds Assumption）¶

実装¶

サンプルデータ¶

statsmodelsによる実装（最尤推定）¶

PyMCによる実装（MCMC）¶

参考¶