順序回帰

[1]:
import numpy as np
import pandas as pd
import scipy.stats as stats

from statsmodels.miscmodels.ordinal_model import OrderedModel

UCLAのウェブサイトからstataデータファイルを読み込みます。このノートブックは、https://stats.idre.ucla.edu/r/dae/ordinal-logistic-regression/(UCLAのRノートブック)に触発されています。

[2]:
url = "https://stats.idre.ucla.edu/stat/data/ologit.dta"
data_student = pd.read_stata(url)
[3]:
data_student.head(5)
[3]:
志願 片親 公立 GPA
0 非常に可能性が高い 0 0 3.26
1 やや可能性が高い 1 0 3.21
2 可能性が低い 1 1 3.94
3 やや可能性が高い 0 0 2.81
4 やや可能性が高い 0 0 2.53
[4]:
data_student.dtypes
[4]:
apply     category
pared         int8
public        int8
gpa        float32
dtype: object
[5]:
data_student['apply'].dtype
[5]:
CategoricalDtype(categories=['unlikely', 'somewhat likely', 'very likely'], ordered=True, categories_dtype=object)

このデータセットは、3つの外生変数(- 0から4までの浮動小数点数であるGPA(gpa)、- 片親が大学院に進学したかどうかを示す2値変数であるpared、- 学生の現在の学部が公立か私立かを示す2値変数であるpublic)を考慮した際に、学部生が大学院に進学する確率に関するものです。

apply(目的変数)は、順序付きカテゴリを持つカテゴリ変数です。unlikely < somewhat likely < very likely。これはカテゴリ型pd.Serieです。NumPy配列よりも好ましいです。

このモデルは、観測できないが外生変数から計算できる潜在的な数値変数\(y_{latent}\)に基づいています。さらに、この\(y_{latent}\)を使用して、観測可能な\(y\)を定義できます。

詳細については、OrderedModelのドキュメント、UCLAのウェブページ、またはこの書籍を参照してください。

プロビット順序回帰:

[6]:
mod_prob = OrderedModel(data_student['apply'],
                        data_student[['pared', 'public', 'gpa']],
                        distr='probit')

res_prob = mod_prob.fit(method='bfgs')
res_prob.summary()
Optimization terminated successfully.
         Current function value: 0.896869
         Iterations: 17
         Function evaluations: 21
         Gradient evaluations: 21
[6]:
OrderedModelの結果
従属変数 志願 対数尤度 -358.75
モデル OrderedModel AIC 727.5
方法 最尤法 BIC 747.5
日付 2024年10月3日(木)
時刻 15:51:00
観測数 400
残差の自由度 395
モデルの自由度 3
係数 標準誤差 z P>|z| [0.025 0.975]
片親 0.5981 0.158 3.789 0.000 0.289 0.908
公立 0.0102 0.173 0.059 0.953 -0.329 0.349
GPA 0.3582 0.157 2.285 0.022 0.051 0.665
可能性が低い/やや可能性が高い 1.2968 0.468 2.774 0.006 0.381 2.213
やや可能性が高い/非常に可能性が高い 0.1873 0.074 2.530 0.011 0.042 0.332

私たちのモデルでは、3つの外生変数(ドキュメントの表記に従うなら\(\beta\))があるので、推定する必要がある係数は3つです。

これらの3つの推定値とその標準誤差は、要約表から取得できます。

目的変数には3つのカテゴリ(unlikelysomewhat likelyvery likely)があるので、推定する必要がある閾値は2つあります。OrderedModel.transform_threshold_paramsメソッドのドキュメントで説明されているように、最初に推定される閾値は実際の値であり、他のすべての閾値は累積指数関数的増分の観点から表されます。実際の閾値は、次のように計算できます。

[7]:
num_of_thresholds = 2
mod_prob.transform_threshold_params(res_prob.params[-num_of_thresholds:])
[7]:
array([      -inf, 1.29684541, 2.50285885,        inf])

ロジット順序回帰:

[8]:
mod_log = OrderedModel(data_student['apply'],
                        data_student[['pared', 'public', 'gpa']],
                        distr='logit')

res_log = mod_log.fit(method='bfgs', disp=False)
res_log.summary()
[8]:
OrderedModelの結果
従属変数 志願 対数尤度 -358.51
モデル OrderedModel AIC 727.0
方法 最尤法 BIC 747.0
日付 2024年10月3日(木)
時刻 15:51:00
観測数 400
残差の自由度 395
モデルの自由度 3
係数 標準誤差 z P>|z| [0.025 0.975]
片親 1.0476 0.266 3.942 0.000 0.527 1.569
公立 -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
可能性が低い/やや可能性が高い 2.2035 0.780 2.827 0.005 0.676 3.731
やや可能性が高い/非常に可能性が高い 0.7398 0.080 9.236 0.000 0.583 0.897
[9]:
predicted = res_log.model.predict(res_log.params, exog=data_student[['pared', 'public', 'gpa']])
predicted
[9]:
array([[0.54884071, 0.35932276, 0.09183653],
       [0.30558191, 0.47594216, 0.21847593],
       [0.22938356, 0.47819057, 0.29242587],
       ...,
       [0.69380357, 0.25470075, 0.05149568],
       [0.54884071, 0.35932276, 0.09183653],
       [0.50896793, 0.38494062, 0.10609145]])
[10]:
pred_choice = predicted.argmax(1)
print('Fraction of correct choice predictions')
print((np.asarray(data_student['apply'].values.codes) == pred_choice).mean())
Fraction of correct choice predictions
0.5775

カスタム累積cLogLog分布による順序回帰:

logit回帰とprobit回帰に加えて、SciPy.statsパッケージの任意の連続分布をdistr引数に使用できます。あるいは、rv_continuousからサブクラスを作成し、いくつかのメソッドを実装するだけで、独自の分布を定義することもできます。

[11]:
# using a SciPy distribution
res_exp = OrderedModel(data_student['apply'],
                           data_student[['pared', 'public', 'gpa']],
                           distr=stats.expon).fit(method='bfgs', disp=False)
res_exp.summary()
[11]:
OrderedModelの結果
従属変数 志願 対数尤度 -360.84
モデル OrderedModel AIC 731.7
方法 最尤法 BIC 751.6
日付 2024年10月3日(木)
時刻 15:51:01
観測数 400
残差の自由度 395
モデルの自由度 3
係数 標準誤差 z P>|z| [0.025 0.975]
片親 0.4690 0.117 4.021 0.000 0.240 0.698
公立 -0.1308 0.149 -0.879 0.379 -0.422 0.161
GPA 0.2198 0.134 1.638 0.101 -0.043 0.483
可能性が低い/やや可能性が高い 1.5370 0.405 3.792 0.000 0.742 2.332
やや可能性が高い/非常に可能性が高い 0.4082 0.093 4.403 0.000 0.226 0.590
[12]:
# minimal definition of a custom scipy distribution.
class CLogLog(stats.rv_continuous):
    def _ppf(self, q):
        return np.log(-np.log(1 - q))

    def _cdf(self, x):
        return 1 - np.exp(-np.exp(x))


cloglog = CLogLog()

# definition of the model and fitting
res_cloglog = OrderedModel(data_student['apply'],
                           data_student[['pared', 'public', 'gpa']],
                           distr=cloglog).fit(method='bfgs', disp=False)
res_cloglog.summary()
[12]:
OrderedModelの結果
従属変数 志願 対数尤度 -359.75
モデル OrderedModel AIC 729.5
方法 最尤法 BIC 749.5
日付 2024年10月3日(木)
時刻 15:51:01
観測数 400
残差の自由度 395
モデルの自由度 3
係数 標準誤差 z P>|z| [0.025 0.975]
片親 0.5167 0.161 3.202 0.001 0.200 0.833
公立 0.1081 0.168 0.643 0.520 -0.221 0.438
GPA 0.3344 0.154 2.168 0.030 0.032 0.637
可能性が低い/やや可能性が高い 0.8705 0.455 1.912 0.056 -0.022 1.763
やや可能性が高い/非常に可能性が高い 0.0989 0.071 1.384 0.167 -0.041 0.239

数式を使用 - 内生変数の処理

Pandasの順序付きカテゴリ変数と数値変数は、数式における従属変数としてサポートされています。他の型はValueErrorを発生させます。

[13]:
modf_logit = OrderedModel.from_formula("apply ~ 0 + pared + public + gpa", data_student,
                                      distr='logit')
resf_logit = modf_logit.fit(method='bfgs')
resf_logit.summary()
Optimization terminated successfully.
         Current function value: 0.896281
         Iterations: 22
         Function evaluations: 24
         Gradient evaluations: 24
[13]:
OrderedModelの結果
従属変数 志願 対数尤度 -358.51
モデル OrderedModel AIC 727.0
方法 最尤法 BIC 747.0
日付 2024年10月3日(木)
時刻 15:51:02
観測数 400
残差の自由度 395
モデルの自由度 3
係数 標準誤差 z P>|z| [0.025 0.975]
片親 1.0476 0.266 3.942 0.000 0.527 1.569
公立 -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
可能性が低い/やや可能性が高い 2.2035 0.780 2.827 0.005 0.676 3.731
やや可能性が高い/非常に可能性が高い 0.7398 0.080 9.236 0.000 0.583 0.897

従属変数に数値コードを使用することはサポートされていますが、カテゴリレベルの名前が失われます。レベルと名前は、数式を使用しない場合と同様に、アルファベット順にソートされた従属変数の固有値に対応します。

[14]:
data_student["apply_codes"] = data_student['apply'].cat.codes * 2 + 5
data_student["apply_codes"].head()
[14]:
0    9
1    7
2    5
3    7
4    7
Name: apply_codes, dtype: int8
[15]:
OrderedModel.from_formula("apply_codes ~ 0 + pared + public + gpa", data_student,
                          distr='logit').fit().summary()
Optimization terminated successfully.
         Current function value: 0.896281
         Iterations: 421
         Function evaluations: 663
[15]:
OrderedModelの結果
従属変数 apply_codes 対数尤度 -358.51
モデル OrderedModel AIC 727.0
方法 最尤法 BIC 747.0
日付 2024年10月3日(木)
時刻 15:51:02
観測数 400
残差の自由度 395
モデルの自由度 3
係数 標準誤差 z P>|z| [0.025 0.975]
片親 1.0477 0.266 3.942 0.000 0.527 1.569
公立 -0.0587 0.298 -0.197 0.844 -0.642 0.525
GPA 0.6157 0.261 2.362 0.018 0.105 1.127
5.0/7.0 2.2033 0.780 2.826 0.005 0.675 3.731
7.0/9.0 0.7398 0.080 9.236 0.000 0.583 0.897
[16]:
resf_logit.predict(data_student.iloc[:5])
[16]:
0 1 2
0 0.548841 0.359323 0.091837
1 0.305582 0.475942 0.218476
2 0.229384 0.478191 0.292426
3 0.616118 0.312690 0.071191
4 0.656003 0.283398 0.060599

従属変数として文字列値を直接使用すると、ValueErrorが発生します。

[17]:
data_student["apply_str"] = np.asarray(data_student["apply"])
data_student["apply_str"].head()
[17]:
0        very likely
1    somewhat likely
2           unlikely
3    somewhat likely
4    somewhat likely
Name: apply_str, dtype: object
[18]:
data_student.apply_str = pd.Categorical(data_student.apply_str, ordered=True)
data_student.public = data_student.public.astype(float)
data_student.pared = data_student.pared.astype(float)
[19]:
OrderedModel.from_formula("apply_str ~ 0 + pared + public + gpa", data_student,
                          distr='logit')
[19]:
<statsmodels.miscmodels.ordinal_model.OrderedModel at 0x7fecb8fbd420>

数式を使用 - モデルに定数項を含めない

OrderedModelのパラメータ化では、モデルに**定数項**を含めない必要があります。明示的にも暗示的にも含めないでください。定数項はすべての閾値をシフトすることと等価であるため、個別に識別されません。

説明変数の中にカテゴリ変数(またはスプライン)がある場合、Patsyの数式仕様では、明示的または暗示的な定数項のない計画行列は許可されません。回避策として、statsmodelsは明示的な切片を除去します。

したがって、切片のない計画行列を取得するための2つの有効なケースがあります。

  • モデルに明示的および暗示的な切片を含めないように指定します。これは、モデルに数値変数のみが含まれている場合に可能です。

  • 明示的な切片を含むモデルを指定します。statsmodelsはこれを除去します。

暗示的な切片を持つモデルは過剰パラメータ化され、パラメータ推定値は完全に識別されず、cov_paramsは反転せず、標準誤差にはNaNが含まれる可能性があります。

次に、追加のカテゴリ変数を含む例を見てみましょう。

[20]:
nobs = len(data_student)
data_student["dummy"] = (np.arange(nobs) < (nobs / 2)).astype(float)

**明示的な切片**(除去されます)

「1 +」は、Patsyのデフォルトであるため、ここでは冗長です。

[21]:
modfd_logit = OrderedModel.from_formula("apply ~ 1 + pared + public + gpa + C(dummy)", data_student,
                                      distr='logit')
resfd_logit = modfd_logit.fit(method='bfgs')
print(resfd_logit.summary())
Optimization terminated successfully.
         Current function value: 0.896247
         Iterations: 26
         Function evaluations: 28
         Gradient evaluations: 28
                             OrderedModel Results
==============================================================================
Dep. Variable:                  apply   Log-Likelihood:                -358.50
Model:                   OrderedModel   AIC:                             729.0
Method:            Maximum Likelihood   BIC:                             752.9
Date:                Thu, 03 Oct 2024
Time:                        15:51:03
No. Observations:                 400
Df Residuals:                     394
Df Model:                           4
===============================================================================================
                                  coef    std err          z      P>|z|      [0.025      0.975]
-----------------------------------------------------------------------------------------------
C(dummy)[T.1.0]                 0.0326      0.198      0.164      0.869      -0.356       0.421
pared                           1.0489      0.266      3.945      0.000       0.528       1.570
public                         -0.0589      0.298     -0.198      0.843      -0.643       0.525
gpa                             0.6153      0.261      2.360      0.018       0.104       1.126
unlikely/somewhat likely        2.2183      0.785      2.826      0.005       0.680       3.757
somewhat likely/very likely     0.7398      0.080      9.237      0.000       0.583       0.897
===============================================================================================
[22]:
modfd_logit.k_vars
[22]:
4
[23]:
modfd_logit.k_constant
[23]:
0

**暗示的な切片**は、過剰パラメータ化されたモデルを作成します

数式で「0 +」を指定すると、明示的な切片が削除されます。ただし、カテゴリエンコーディングが変更され、暗示的な切片が含まれるようになります。この例では、作成されたダミー変数C(dummy)[0.0]C(dummy)[1.0]の合計は1になります。

OrderedModel.from_formula("apply ~ 0 + pared + public + gpa + C(dummy)", data_student, distr='logit')

過剰パラメータ化された場合に何が起こるかを確認するために、モデルで定数チェックを明示的に指定することで回避できます。モデルに暗示的な定数がある場合でも、hasconst=Falseを使用します。

2つのダミー変数列のパラメータと最初の閾値は個別に識別されません。これらのパラメータの推定値と標準誤差の可用性は任意であり、環境間で異なる数値の詳細によって異なります。

対数尤度値などのいくつかの要約指標は、収束許容範囲と数値精度内で影響を受けません。予測も可能ですが、推論は利用できません、または有効ではありません。

[24]:
modfd2_logit = OrderedModel.from_formula("apply ~ 0 + pared + public + gpa + C(dummy)", data_student,
                                         distr='logit', hasconst=False)
resfd2_logit = modfd2_logit.fit(method='bfgs')
print(resfd2_logit.summary())
Optimization terminated successfully.
         Current function value: 0.896247
         Iterations: 24
         Function evaluations: 26
         Gradient evaluations: 26
                             OrderedModel Results
==============================================================================
Dep. Variable:                  apply   Log-Likelihood:                -358.50
Model:                   OrderedModel   AIC:                             731.0
Method:            Maximum Likelihood   BIC:                             758.9
Date:                Thu, 03 Oct 2024
Time:                        15:51:03
No. Observations:                 400
Df Residuals:                     393
Df Model:                           5
===============================================================================================
                                  coef    std err          z      P>|z|      [0.025      0.975]
-----------------------------------------------------------------------------------------------
C(dummy)[0.0]                  -0.6834        nan        nan        nan         nan         nan
C(dummy)[1.0]                  -0.6508        nan        nan        nan         nan         nan
pared                           1.0489        nan        nan        nan         nan         nan
public                         -0.0588        nan        nan        nan         nan         nan
gpa                             0.6153        nan        nan        nan         nan         nan
unlikely/somewhat likely        1.5349        nan        nan        nan         nan         nan
somewhat likely/very likely     0.7398        nan        nan        nan         nan         nan
===============================================================================================
/opt/hostedtoolcache/Python/3.10.15/x64/lib/python3.10/site-packages/statsmodels/base/model.py:595: HessianInversionWarning: Inverting hessian failed, no bse or cov_params available
  warnings.warn('Inverting hessian failed, no bse or cov_params '
[25]:
resfd2_logit.predict(data_student.iloc[:5])
[25]:
0 1 2
0 0.544858 0.361972 0.093170
1 0.301918 0.476667 0.221416
2 0.226434 0.477700 0.295867
3 0.612254 0.315481 0.072264
4 0.652280 0.286188 0.061532
[26]:
resf_logit.predict()
[26]:
array([[0.54884071, 0.35932276, 0.09183653],
       [0.30558191, 0.47594216, 0.21847593],
       [0.22938356, 0.47819057, 0.29242587],
       ...,
       [0.69380357, 0.25470075, 0.05149568],
       [0.54884071, 0.35932276, 0.09183653],
       [0.50896793, 0.38494062, 0.10609145]])

二項モデルとロジットの比較

従属順序カテゴリ変数に2つのレベルしかない場合、モデルはロジットモデルでも推定できます。

この場合、定数のパラメータ化を除いて、モデルは(理論的には)同一です。ロジットは他のほとんどのモデルと同様に、一般的に切片を必要とします。これはOrderedModelの閾値パラメータに対応しますが、符号が反対です。

実装は異なり、同じ結果統計と推定後の機能のすべてが利用できるわけではありません。推定されたパラメータと他の結果統計は、主に推定の収束許容範囲に基づいて異なります。

[27]:
from statsmodels.discrete.discrete_model import Logit
from statsmodels.tools.tools import add_constant

データから中央のカテゴリを削除し、2つの極端なカテゴリを保持します。

[28]:
mask_drop = data_student['apply'] == "somewhat likely"
data2 = data_student.loc[~mask_drop, :].copy()
# we need to remove the category also from the Categorical Index
data2['apply'] = data2['apply'].cat.remove_categories("somewhat likely")
data2["apply"].head()
[28]:
0     very likely
2        unlikely
5        unlikely
8        unlikely
10       unlikely
Name: apply, dtype: category
Categories (2, object): ['unlikely' < 'very likely']
[29]:
mod_log = OrderedModel(data2['apply'],
                        data2[['pared', 'public', 'gpa']],
                        distr='logit')

res_log = mod_log.fit(method='bfgs', disp=False)
res_log.summary()
[29]:
OrderedModelの結果
従属変数 志願 対数尤度 -102.87
モデル OrderedModel AIC 213.7
方法 最尤法 BIC 228.0
日付 2024年10月3日(木)
時刻 15:51:04
観測数 260
残差の自由度 256
モデルの自由度 3
係数 標準誤差 z P>|z| [0.025 0.975]
片親 1.2861 0.438 2.934 0.003 0.427 2.145
公立 0.4014 0.444 0.903 0.366 -0.470 1.272
GPA 0.7854 0.489 1.605 0.108 -0.174 1.744
可能性が低い/非常に可能性が高い 4.4147 1.485 2.974 0.003 1.505 7.324

ロジットモデルにはデフォルトで定数項がありません。説明変数に追加する必要があります。

結果は、主に推定における収束許容範囲に起因する数値精度まで、ロジットと順序モデルの間で本質的に同一です。

唯一の違いは定数の符号です。ロジットとOrderedModelでは、定数の符号が反対です。これは、計画行列に定数列を含めるのではなく、OrderedModelでカットポイントの観点からパラメータ化されるためです。

[30]:
ex = add_constant(data2[['pared', 'public', 'gpa']], prepend=False)
mod_logit = Logit(data2['apply'].cat.codes, ex)

res_logit = mod_logit.fit(method='bfgs', disp=False)
[31]:
res_logit.summary()
[31]:
ロジット回帰の結果
従属変数 y 観測数 260
モデル ロジット 残差の自由度 256
方法 最尤推定 モデルの自由度 3
日付 2024年10月3日(木) 擬似R二乗 0.07842
時刻 15:51:04 対数尤度 -102.87
収束 True LL-Null -111.62
共分散の種類 頑健でない LLR p値 0.0005560
係数 標準誤差 z P>|z| [0.025 0.975]
片親 1.2861 0.438 2.934 0.003 0.427 2.145
公立 0.4014 0.444 0.903 0.366 -0.470 1.272
GPA 0.7854 0.489 1.605 0.108 -0.174 1.744
定数 -4.4148 1.485 -2.974 0.003 -7.324 -1.505

ロジットと同様に、OrderedModelでも頑健な標準誤差を利用できます。例として、クロスセクションデータがあり、自己相関が適切ではない場合でも、HAC共分散タイプを指定します。

[32]:
res_logit_hac = mod_logit.fit(method='bfgs', disp=False, cov_type="hac", cov_kwds={"maxlags": 2})
res_log_hac = mod_log.fit(method='bfgs', disp=False, cov_type="hac", cov_kwds={"maxlags": 2})
[33]:
res_logit_hac.bse.values - res_log_hac.bse
[33]:
pared                   9.022236e-08
public                 -7.249837e-08
gpa                     7.653233e-08
unlikely/very likely    2.800466e-07
dtype: float64

最終更新日: 2024年10月3日