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

累積ロジットモデル（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) \]

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

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

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

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

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

各カテゴリに属する確率\(P(y = k \mid \mathbf{x})\)は、累積確率の差で表現する。

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

つまり

\[ P(y = k \mid \mathbf{x}) = \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]} \]

と表される。

端のカテゴリの扱い#

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

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

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

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

こうすることで、

  • 最小カテゴリ： \(P(y=1)=P(y \le 1)\)

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

と書ける。

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

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

と表すこともできる

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

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

  • 切片（\(\tau_k\)）は \(k\) ごとに異なる

  • 傾き（\(\beta\)）はすべて共通

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

そして累積ロジットモデルではカテゴリごとに異なる傾き\(\beta_k\)を持つとされる…が、共通の傾きを持つように書く文献のほうが多い

実装#

サンプルデータ#

  • 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)
apply pared public gpa
0 very likely 0 0 3.26
1 somewhat likely 1 0 3.21
2 unlikely 1 1 3.94
3 somewhat likely 0 0 2.81
4 somewhat likely 0 0 2.53

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
OrderedModel Results
Dep. Variable: apply Log-Likelihood: -358.51
Model: OrderedModel AIC: 727.0
Method: Maximum Likelihood BIC: 747.0
Date: Sat, 14 Feb 2026
Time: 09:40:06
No. Observations: 400
Df Residuals: 395
Df Model: 3
coef std err z P>|z| [0.025 0.975]
pared 1.0476 0.266 3.942 0.000 0.527 1.569
public -0.0586 0.298 -0.197 0.844 -0.642 0.525
gpa 0.6158 0.261 2.363 0.018 0.105 1.127
unlikely/somewhat likely 2.2035 0.780 2.827 0.005 0.676 3.731
somewhat likely/very likely 0.7398 0.080 9.236 0.000 0.583 0.897

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]



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),
})













  
    

      
      exog
      EAP
    

  
  
    

      0
      pared
      0.9995
    

    

      1
      public
      -0.0572
    

    

      2
      gpa
      0.5796
    

  













参考#







    
    

                

              

              
              
              
              
                
              
            

            
            
              
                



  

  

     Contents