序數迴歸¶
[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]:
| apply | pared | public | 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)
此資料集關於大學生在給定三個外生變數的情況下，申請研究所的可能性：- 他們的平均成績 ( gpa )，一個介於 0 和 4 之間的浮點數。 - pared，一個二元變數，表示至少有一位家長上過研究所。 - 以及 public，一個二元變數，表示學生目前就讀的大學是公立還是私立。
apply，目標變數是具有排序類別的類別變數：unlikely < somewhat likely < very likely。它是一個 pd.Serie 的類別類型，這比 NumPy 陣列更佳。
此模型基於一個數值潛在變數 \(y_{latent}\)，我們無法觀察到，但可以藉由外生變數計算。此外，我們可以利用這個 \(y_{latent}\) 來定義我們可以觀察到的 \(y\)。
如需更多詳細資訊，請參閱 OrderedModel 的文件、UCLA 網頁 或這本書籍。
Probit 序數迴歸：¶
[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]:
| 依變數 | apply | 對數概似 | -358.75 |
|---|---|---|---|
| 模型 | OrderedModel | AIC | 727.5 |
| 方法 | 最大概似法 | BIC | 747.5 |
| 日期 | 週四, 2024年10月03日 | ||
| 時間 | 15:51:00 | ||
| 觀測值數量 | 400 | ||
| 殘差自由度 | 395 | ||
| 模型自由度 | 3 |
| coef | 標準誤差 | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| pared | 0.5981 | 0.158 | 3.789 | 0.000 | 0.289 | 0.908 |
| public | 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\)s），因此我們需要估計 3 個係數。
這 3 個估計值及其標準誤差可以在摘要表中檢索。
由於目標變數中有 3 個類別（unlikely、somewhat likely、very likely），我們有兩個閾值需要估計。如 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])
Logit 序數迴歸：¶
[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]:
| 依變數 | apply | 對數概似 | -358.51 |
|---|---|---|---|
| 模型 | OrderedModel | AIC | 727.0 |
| 方法 | 最大概似法 | BIC | 747.0 |
| 日期 | 週四, 2024年10月03日 | ||
| 時間 | 15:51:00 | ||
| 觀測值數量 | 400 | ||
| 殘差自由度 | 395 | ||
| 模型自由度 | 3 |
| coef | 標準誤差 | 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 |
| 不太可能/有點可能 | 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]:
| 依變數 | apply | 對數概似 | -360.84 |
|---|---|---|---|
| 模型 | OrderedModel | AIC | 731.7 |
| 方法 | 最大概似法 | BIC | 751.6 |
| 日期 | 週四, 2024年10月03日 | ||
| 時間 | 15:51:01 | ||
| 觀測值數量 | 400 | ||
| 殘差自由度 | 395 | ||
| 模型自由度 | 3 |
| coef | 標準誤差 | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| pared | 0.4690 | 0.117 | 4.021 | 0.000 | 0.240 | 0.698 |
| public | -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]:
| 依變數 | apply | 對數概似 | -359.75 |
|---|---|---|---|
| 模型 | OrderedModel | AIC | 729.5 |
| 方法 | 最大概似法 | BIC | 749.5 |
| 日期 | 週四, 2024年10月03日 | ||
| 時間 | 15:51:01 | ||
| 觀測值數量 | 400 | ||
| 殘差自由度 | 395 | ||
| 模型自由度 | 3 |
| coef | 標準誤差 | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| pared | 0.5167 | 0.161 | 3.202 | 0.001 | 0.200 | 0.833 |
| public | 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]:
| 依變數 | apply | 對數概似 | -358.51 |
|---|---|---|---|
| 模型 | OrderedModel | AIC | 727.0 |
| 方法 | 最大概似法 | BIC | 747.0 |
| 日期 | 週四, 2024年10月03日 | ||
| 時間 | 15:51:02 | ||
| 觀測值數量 | 400 | ||
| 殘差自由度 | 395 | ||
| 模型自由度 | 3 |
| coef | 標準誤差 | 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 |
| 不太可能/有點可能 | 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]:
| 依變數 | apply_codes | 對數概似 | -358.51 |
|---|---|---|---|
| 模型 | OrderedModel | AIC | 727.0 |
| 方法 | 最大概似法 | BIC | 747.0 |
| 日期 | 週四, 2024年10月03日 | ||
| 時間 | 15:51:02 | ||
| 觀測值數量 | 400 | ||
| 殘差自由度 | 395 | ||
| 模型自由度 | 3 |
| coef | 標準誤差 | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| pared | 1.0477 | 0.266 | 3.942 | 0.000 | 0.527 | 1.569 |
| public | -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 會移除一個顯式截距。
因此，有兩種有效情況可以取得沒有截距的設計矩陣。
指定一個沒有顯式和隱式截距的模型，如果模型中只有數值變數，這是可能的。
指定一個帶有顯式截距的模型，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，即使模型具有隱式常數。
兩個虛擬變數欄位的參數和第一個閾值無法單獨識別。這些參數的估計值和標準誤差的可用性是任意的，並且取決於各環境之間不同的數值詳細資訊。
某些摘要度量（如對數概似值）不受此影響，在收斂容差和數值精度的範圍內。應該也可以進行預測。但是，無法進行推論，或者推論無效。
[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]])
二元模型與 Logit 的比較¶
如果依變數的排序類別變數只有兩個層級，則也可以使用 Logit 模型估計該模型。
在這種情況下，這些模型（理論上）是相同的，除了常數的參數化之外。Logit 與大多數其他模型一樣，通常需要截距。這對應於 OrderedModel 中的閾值參數，但符號相反。
實作方式不同，並且並非所有相同的結果統計資料和後估計功能都可用。估計的參數和其他結果統計資料的主要差異基於最佳化的收斂容差。
[27]:
from statsmodels.discrete.discrete_model import Logit
from statsmodels.tools.tools import add_constant
我們從資料中刪除中間類別，並保留兩個極端類別。
[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]:
| 依變數 | apply | 對數概似 | -102.87 |
|---|---|---|---|
| 模型 | OrderedModel | AIC | 213.7 |
| 方法 | 最大概似法 | BIC | 228.0 |
| 日期 | 週四, 2024年10月03日 | ||
| 時間 | 15:51:04 | ||
| 觀測值數量 | 260 | ||
| 殘差自由度 | 256 | ||
| 模型自由度 | 3 |
| coef | 標準誤差 | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| pared | 1.2861 | 0.438 | 2.934 | 0.003 | 0.427 | 2.145 |
| public | 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 |
預設情況下，Logit 模型沒有常數，我們必須將其新增到我們的解釋變數中。
Logit 和排序模型之間的結果基本上是相同的，在數值精度方面主要來自估計中的收斂容差。
唯一的區別在於常數的符號，Logit 和 OrdereModel 的常數符號相反。這是 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 |
|---|---|---|---|
| 模型 | Logit | 殘差自由度 | 256 |
| 方法 | MLE | 模型自由度 | 3 |
| 日期 | 週四, 2024年10月03日 | 偽 R 平方 | 0.07842 |
| 時間 | 15:51:04 | 對數概似 | -102.87 |
| 已收斂 | True | LL-Null | -111.62 |
| 共變異數類型 | 非穩健 | LLR p 值 | 0.0005560 |
| coef | 標準誤差 | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| pared | 1.2861 | 0.438 | 2.934 | 0.003 | 0.427 | 2.145 |
| public | 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 |
| const | -4.4148 | 1.485 | -2.974 | 0.003 | -7.324 | -1.505 |
在 OrderedModel 中，也可以使用與 discrete.Logit 相同的方式使用穩健標準誤差。作為範例，我們指定 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月03日