[1]:
import numpy as np
import pandas as pd
from scipy import stats
from statsmodels.miscmodels.ordinal_model import OrderedModel
从 UCLA 网站加载 Stata 数据文件。本笔记本参考了 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 | 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 |
[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=str)
这个数据集使用三个外生变量，研究本科生申请研究生院的可能性：
-
平均绩点
gpa，取值为 0 到 4 之间的浮点数。 -
二元变量
pared，表示父母中是否至少有一人上过研究生院。 -
二元变量
public，表示学生目前就读的本科学校是公立还是私立。
目标变量 apply 是具有顺序的分类变量，类别关系为 unlikely < somewhat likely < very likely。它使用分类类型的 pandas Series 表示；与 NumPy 数组相比，优先推荐这种表示方式。
模型建立在数值型潜变量 ylatent 上。我们无法直接观测这个变量，但可以通过外生变量计算它，并用它定义可以观测到的 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]:
| Dep. Variable: | apply | Log-Likelihood: | -358.75 |
|---|---|---|---|
| Model: | OrderedModel | AIC: | 727.5 |
| Method: | Maximum Likelihood | BIC: | 747.5 |
| Date: | Thu, 27 Aug 2026 | ||
| Time: | 06:50:50 | ||
| No. Observations: | 400 | ||
| Df Residuals: | 395 | ||
| Df Model: | 3 |
| coef | std err | 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 |
| unlikely/somewhat likely | 1.2968 | 0.468 | 2.774 | 0.006 | 0.381 | 2.213 |
| somewhat likely/very likely | 0.1873 | 0.074 | 2.530 | 0.011 | 0.042 | 0.332 |
这个模型有三个外生变量；沿用文档中的记号，它们对应的系数用 β 表示。因此，需要估计三个系数。
这三个系数的估计值及其标准误，可以在下面的汇总表中找到。
目标变量有三个类别（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]:
| Dep. Variable: | apply | Log-Likelihood: | -358.51 |
|---|---|---|---|
| Model: | OrderedModel | AIC: | 727.0 |
| Method: | Maximum Likelihood | BIC: | 747.0 |
| Date: | Thu, 27 Aug 2026 | ||
| Time: | 06:50:51 | ||
| 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 |
[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.50896794, 0.38494062, 0.10609145]], shape=(400, 3))
[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 回归，distr 参数还可以使用 scipy.stats 中的任意连续分布。也可以通过继承 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]:
| Dep. Variable: | apply | Log-Likelihood: | -360.84 |
|---|---|---|---|
| Model: | OrderedModel | AIC: | 731.7 |
| Method: | Maximum Likelihood | BIC: | 751.6 |
| Date: | Thu, 27 Aug 2026 | ||
| Time: | 06:50:51 | ||
| No. Observations: | 400 | ||
| Df Residuals: | 395 | ||
| Df Model: | 3 |
| coef | std err | 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 |
| unlikely/somewhat likely | 1.5370 | 0.405 | 3.792 | 0.000 | 0.742 | 2.332 |
| somewhat likely/very likely | 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]:
| Dep. Variable: | apply | Log-Likelihood: | -359.75 |
|---|---|---|---|
| Model: | OrderedModel | AIC: | 729.5 |
| Method: | Maximum Likelihood | BIC: | 749.5 |
| Date: | Thu, 27 Aug 2026 | ||
| Time: | 06:50:51 | ||
| No. Observations: | 400 | ||
| Df Residuals: | 395 | ||
| Df Model: | 3 |
| coef | std err | 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 |
| unlikely/somewhat likely | 0.8705 | 0.455 | 1.912 | 0.056 | -0.022 | 1.763 |
| somewhat likely/very likely | 0.0989 | 0.071 | 1.384 | 0.167 | -0.041 | 0.239 |
使用公式：处理 endog
公式支持把 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]:
| Dep. Variable: | apply | Log-Likelihood: | -358.51 |
|---|---|---|---|
| Model: | OrderedModel | AIC: | 727.0 |
| Method: | Maximum Likelihood | BIC: | 747.0 |
| Date: | Thu, 27 Aug 2026 | ||
| Time: | 06:50:51 | ||
| 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 |
也支持用数字编码表示因变量，但这样会丢失类别层级的名称。与不使用公式时一样，层级及其名称对应于因变量的唯一值，并按字母数字顺序排列。
[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]:
| Dep. Variable: | apply_codes | Log-Likelihood: | -358.51 |
|---|---|---|---|
| Model: | OrderedModel | AIC: | 727.0 |
| Method: | Maximum Likelihood | BIC: | 747.0 |
| Date: | Thu, 27 Aug 2026 | ||
| Time: | 06:50:52 | ||
| No. Observations: | 400 | ||
| Df Residuals: | 395 | ||
| Df Model: | 3 |
| coef | std err | 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: str
[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 0x7f2855746360>
使用公式：模型中不能包含常数项
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, 27 Aug 2026
Time: 06:50:52
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, 27 Aug 2026
Time: 06:50:52
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 615.081 -0.001 0.999 -1206.220 1204.853
C(dummy)[1.0] -0.6508 615.081 -0.001 0.999 -1206.187 1204.886
pared 1.0489 0.266 3.944 0.000 0.528 1.570
public -0.0588 0.298 -0.197 0.844 -0.643 0.525
gpa 0.6153 0.261 2.360 0.018 0.104 1.126
unlikely/somewhat likely 1.5349 615.081 0.002 0.998 -1204.003 1207.072
somewhat likely/very likely 0.7398 0.080 9.237 0.000 0.583 0.897
===============================================================================================
[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.50896794, 0.38494062, 0.10609145]], shape=(400, 3))
二元模型与 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, str): ['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]:
| Dep. Variable: | apply | Log-Likelihood: | -102.87 |
|---|---|---|---|
| Model: | OrderedModel | AIC: | 213.7 |
| Method: | Maximum Likelihood | BIC: | 228.0 |
| Date: | Thu, 27 Aug 2026 | ||
| Time: | 06:50:53 | ||
| No. Observations: | 260 | ||
| Df Residuals: | 256 | ||
| Df Model: | 3 |
| coef | std err | 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 |
| unlikely/very likely | 4.4147 | 1.485 | 2.974 | 0.003 | 1.505 | 7.324 |
Logit 模型默认不会自动加入常数项，因此需要把常数项添加到解释变量中。
Logit 与有序模型的结果基本一致；细微的数值差异主要来自估计过程中的收敛容差。
唯一的区别是常数项的符号：Logit 与 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]:
| Dep. Variable: | y | No. Observations: | 260 |
|---|---|---|---|
| Model: | Logit | Df Residuals: | 256 |
| Method: | MLE | Df Model: | 3 |
| Date: | Thu, 27 Aug 2026 | Pseudo R-squ.: | 0.07842 |
| Time: | 06:50:53 | Log-Likelihood: | -102.87 |
| converged: | True | LL-Null: | -111.62 |
| Covariance Type: | nonrobust | LLR p-value: | 0.0005560 |
| coef | std err | 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
© 版权声明
文章版权归作者所有，未经允许请勿转载。
THE END
暂无评论内容