多项式Logit模型Python和Stata不同的结果
Multinomial Logit model Python and Stata different results
我正在尝试使用 python 和 stata 构建多项式 logit 模型。我的数据如下：
ses_type prog_type read write math prog ses
0 low Diploma 39.2 40.2 46.2 0 0
1 middle general 39.2 38.2 46.2 1 1
2 high Diploma 44.5 44.5 49.5 0 2
3 low Diploma 43.0 43.0 48.0 0 0
4 middle Diploma 44.5 36.5 45.5 0 1
5 high general 47.3 41.3 47.3 1 2
我正在尝试使用 ses 读写和数学 来预测 prog。其中 ses 代表社会经济地位并且是名义变量因此我使用以下命令在 stata 中创建了我的模型：
mlogit prog i.ses read write math, base(2)
Stata输出如下：
Iteration 0: log likelihood = -204.09667
Iteration 1: log likelihood = -171.90258
Iteration 2: log likelihood = -170.13513
Iteration 3: log likelihood = -170.11071
Iteration 4: log likelihood = -170.1107
Multinomial logistic regression Number of obs = 200
LR chi2(10) = 67.97
Prob > chi2 = 0.0000
Log likelihood = -170.1107 Pseudo R2 = 0.1665
------------------------------------------------------------------------------
prog | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
0 |
ses |
1 | .6197969 .5059335 1.23 0.221 -.3718146 1.611408
2 | -.5131952 .6280601 -0.82 0.414 -1.74417 .7177799
|
read | -.0405302 .0289314 -1.40 0.161 -.0972346 .0161742
write | -.0459711 .0270153 -1.70 0.089 -.09892 .0069779
math | -.0990497 .0331576 -2.99 0.003 -.1640373 -.0340621
_cons | 9.544131 1.738404 5.49 0.000 6.136921 12.95134
-------------+----------------------------------------------------------------
1 |
ses |
1 | -.3350861 .4607246 -0.73 0.467 -1.23809 .5679176
2 | -.8687013 .5363968 -1.62 0.105 -1.92002 .182617
|
read | -.0226249 .0264534 -0.86 0.392 -.0744726 .0292228
write | -.011618 .0266782 -0.44 0.663 -.0639063 .0406703
math | -.0591301 .0299996 -1.97 0.049 -.1179283 -.000332
_cons | 5.041193 1.524174 3.31 0.001 2.053866 8.028519
-------------+----------------------------------------------------------------
2 | (base outcome)
------------------------------------------------------------------------------
我尝试使用 python 中的 scikit 学习模块复制相同的结果。以下是代码：
data = pd.read_csv("C://Users/Furqan/Desktop/random_data.csv")
train_x = np.array(data[['read', 'write', 'math','ses ']])
train_y = np.array(data['prog'])
mul_lr = linear_model.LogisticRegression(multi_class='multinomial',
solver='newton-cg').fit(train_x, train_y)
print(mul_lr.intercept_)
print(mul_lr.coef_)
输出值（截距和系数）如下：
[ 4.76438772 0.19347405 -4.95786177]
[[-0.01735513 -0.02731273 -0.04463257 0.01721334]
[-0.00319366 0.00783135 -0.00689664 -0.24480926]
[ 0.02054879 0.01948137 0.05152921 0.22759592]]
原来值不一样。
我的第一个问题是为什么结果往往不同？
我的第二个问题是，在名义预测变量的情况下，我们如何指示 python ses 是一个 指标变量?
编辑：
Link 到数据文件
有几个问题导致 Stata 和 sklearn 结果不同：
- Stata 和 sklearn 中不同的实际预测器
- 拟合参数的不同表示
- 拟合模型时的不同目标函数
我们需要更改所有三个条件才能获得相似的输出。
1。制作虚拟变量
Stata用于线性部分的公式是
prediction = a0 + a1 * [ses==1] + a2 * [ses==2] + a3 * read + a4 * write + a5 * math
Sklearn，反过来，对 ses 的分类性质一无所知，并尝试使用
prediction = a0 + a1 * ses + a3 * read + a4 * write + a5 * math
要启用分类预测，您需要预处理 数据。这是将分类变量包含到 sklearn 逻辑回归中的唯一可能方法。我发现 pd.get_dummies() 是最方便的方法。
以下代码为 ses 创建虚拟变量，然后删除 "low" 级别，这显然对应于示例中的 ses=0：
import pandas as pd, numpy as np
from sklearn import linear_model
data = pd.read_csv("d1.csv", sep='\t')
data.columns = data.columns.str.strip()
raw_x = data.drop('prog', axis=1)
# making the dummies
train_x = pd.get_dummies(raw_x, columns=['ses']).drop('ses_low ', axis=1)
print(train_x.columns)
train_y = data['prog']
mul_lr = linear_model.LogisticRegression(multi_class='multinomial',
solver='newton-cg').fit(train_x, train_y)
reorder = [4, 3, 0, 1, 2] # the order in which coefficents show up in Stata
print(mul_lr.intercept_)
print(mul_lr.coef_[:, reorder])
输出
['read', 'write', 'math', 'ses_high ', 'ses_middle ']
[ 4.67331919 0.19082335 -4.86414254]
[[ 0.47140512 -0.08236331 -0.01909793 -0.02680609 -0.04587383]
[-0.36381476 -0.33294749 -0.0021255 0.00765828 -0.00703075]
[-0.10759035 0.4153108 0.02122343 0.01914781 0.05290458]]
您看到 Python 已成功将 sess 编码为 'ses_high ' 和 'ses_middle '，但未能产生预期的系数。
顺便说一句，我更改了输出中 coef_ 列的顺序，使其看起来像在 Stata 中。
2。重新排列结果
这是因为 Stata 将结果的第三类 (prog=='honors ') 视为 基本结果 ，并从其余参数中减去其所有参数。在 Python 中，您可以通过 运行
重现此内容
print(mul_lr.intercept_ - mul_lr.intercept_[-1])
print((mul_lr.coef_ - mul_lr.coef_[-1])[:, reorder])
这给了你
[9.53746174 5.0549659 0. ]
[[ 0.57899547 -0.4976741 -0.04032136 -0.0459539 -0.09877841]
[-0.25622441 -0.74825829 -0.02334893 -0.01148954 -0.05993533]
[ 0. 0. 0. 0. 0. ]]
您现在可以看到参数现在接近 Stata 给出的参数：
- 截取 Python 中的 (9.53, 5.05) 与 Stata 中的 (9.54, 5.04)
- first-outcome 系数 (0.57, -0.49, ...) 对比 (0.61, -0.51, ...)
- second-outcome 系数 (-0.25, -0.74, ...) 对比 (-0.33, -0.86, ...)
你能看出规律吗？在 sklearn 中，斜率系数比在 Stata 中更小（接近于零）。这不是意外！
3。处理正则化
发生这种情况是因为 sklearn 有意 将斜率系数缩小到 0，方法是在系数上对其最大化的似然函数添加二次惩罚。这使得估计有偏差但更稳定，即使在严重的多重共线性的情况下也是如此。用贝叶斯术语来说，这种正则化对应于所有系数的 zero-mean 高斯先验。您可以了解有关正则化的更多信息 in the wiki.
在sklearn中，这个二次惩罚由正C参数控制：它越小，你得到的正则化程度越高。您可以将其视为每个斜率系数的先验方差。默认值是 C=1，但是你可以把它变大，比如 C=1000000，这意味着几乎没有正则化。在这种情况下，输出几乎与 Stata:
相同
mul_lr2 = linear_model.LogisticRegression(
multi_class='multinomial', solver='newton-cg', C=1000000
).fit(train_x, train_y)
print(mul_lr2.intercept_ - mul_lr2.intercept_[-1])
print((mul_lr2.coef_ - mul_lr2.coef_[-1])[:, reorder])
这给了你
[9.54412644 5.04126452 0. ]
[[ 0.61978951 -0.51320481 -0.04053013 -0.0459711 -0.09904948]
[-0.33508605 -0.86869799 -0.02262518 -0.01161839 -0.05913068]
[ 0. 0. 0. 0. 0. ]]
结果仍然略有不同（如小数点后 5 位），但正则化更少时，差异填充会进一步缩小。
我正在尝试使用 python 和 stata 构建多项式 logit 模型。我的数据如下：
ses_type prog_type read write math prog ses
0 low Diploma 39.2 40.2 46.2 0 0
1 middle general 39.2 38.2 46.2 1 1
2 high Diploma 44.5 44.5 49.5 0 2
3 low Diploma 43.0 43.0 48.0 0 0
4 middle Diploma 44.5 36.5 45.5 0 1
5 high general 47.3 41.3 47.3 1 2
我正在尝试使用 ses 读写和数学 来预测 prog。其中 ses 代表社会经济地位并且是名义变量因此我使用以下命令在 stata 中创建了我的模型：
mlogit prog i.ses read write math, base(2)
Stata输出如下：
Iteration 0: log likelihood = -204.09667
Iteration 1: log likelihood = -171.90258
Iteration 2: log likelihood = -170.13513
Iteration 3: log likelihood = -170.11071
Iteration 4: log likelihood = -170.1107
Multinomial logistic regression Number of obs = 200
LR chi2(10) = 67.97
Prob > chi2 = 0.0000
Log likelihood = -170.1107 Pseudo R2 = 0.1665
------------------------------------------------------------------------------
prog | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
0 |
ses |
1 | .6197969 .5059335 1.23 0.221 -.3718146 1.611408
2 | -.5131952 .6280601 -0.82 0.414 -1.74417 .7177799
|
read | -.0405302 .0289314 -1.40 0.161 -.0972346 .0161742
write | -.0459711 .0270153 -1.70 0.089 -.09892 .0069779
math | -.0990497 .0331576 -2.99 0.003 -.1640373 -.0340621
_cons | 9.544131 1.738404 5.49 0.000 6.136921 12.95134
-------------+----------------------------------------------------------------
1 |
ses |
1 | -.3350861 .4607246 -0.73 0.467 -1.23809 .5679176
2 | -.8687013 .5363968 -1.62 0.105 -1.92002 .182617
|
read | -.0226249 .0264534 -0.86 0.392 -.0744726 .0292228
write | -.011618 .0266782 -0.44 0.663 -.0639063 .0406703
math | -.0591301 .0299996 -1.97 0.049 -.1179283 -.000332
_cons | 5.041193 1.524174 3.31 0.001 2.053866 8.028519
-------------+----------------------------------------------------------------
2 | (base outcome)
------------------------------------------------------------------------------
我尝试使用 python 中的 scikit 学习模块复制相同的结果。以下是代码：
data = pd.read_csv("C://Users/Furqan/Desktop/random_data.csv")
train_x = np.array(data[['read', 'write', 'math','ses ']])
train_y = np.array(data['prog'])
mul_lr = linear_model.LogisticRegression(multi_class='multinomial',
solver='newton-cg').fit(train_x, train_y)
print(mul_lr.intercept_)
print(mul_lr.coef_)
输出值（截距和系数）如下：
[ 4.76438772 0.19347405 -4.95786177]
[[-0.01735513 -0.02731273 -0.04463257 0.01721334]
[-0.00319366 0.00783135 -0.00689664 -0.24480926]
[ 0.02054879 0.01948137 0.05152921 0.22759592]]
原来值不一样。
我的第一个问题是为什么结果往往不同？
我的第二个问题是，在名义预测变量的情况下，我们如何指示 python ses 是一个 指标变量?
编辑：
Link 到数据文件
有几个问题导致 Stata 和 sklearn 结果不同：
- Stata 和 sklearn 中不同的实际预测器
- 拟合参数的不同表示
- 拟合模型时的不同目标函数
我们需要更改所有三个条件才能获得相似的输出。
1。制作虚拟变量
Stata用于线性部分的公式是
prediction = a0 + a1 * [ses==1] + a2 * [ses==2] + a3 * read + a4 * write + a5 * math
Sklearn，反过来，对 ses 的分类性质一无所知，并尝试使用
prediction = a0 + a1 * ses + a3 * read + a4 * write + a5 * math
要启用分类预测，您需要预处理 数据。这是将分类变量包含到 sklearn 逻辑回归中的唯一可能方法。我发现 pd.get_dummies() 是最方便的方法。
以下代码为 ses 创建虚拟变量，然后删除 "low" 级别，这显然对应于示例中的 ses=0：
import pandas as pd, numpy as np
from sklearn import linear_model
data = pd.read_csv("d1.csv", sep='\t')
data.columns = data.columns.str.strip()
raw_x = data.drop('prog', axis=1)
# making the dummies
train_x = pd.get_dummies(raw_x, columns=['ses']).drop('ses_low ', axis=1)
print(train_x.columns)
train_y = data['prog']
mul_lr = linear_model.LogisticRegression(multi_class='multinomial',
solver='newton-cg').fit(train_x, train_y)
reorder = [4, 3, 0, 1, 2] # the order in which coefficents show up in Stata
print(mul_lr.intercept_)
print(mul_lr.coef_[:, reorder])
输出
['read', 'write', 'math', 'ses_high ', 'ses_middle ']
[ 4.67331919 0.19082335 -4.86414254]
[[ 0.47140512 -0.08236331 -0.01909793 -0.02680609 -0.04587383]
[-0.36381476 -0.33294749 -0.0021255 0.00765828 -0.00703075]
[-0.10759035 0.4153108 0.02122343 0.01914781 0.05290458]]
您看到 Python 已成功将 sess 编码为 'ses_high ' 和 'ses_middle '，但未能产生预期的系数。
顺便说一句，我更改了输出中 coef_ 列的顺序，使其看起来像在 Stata 中。
2。重新排列结果
这是因为 Stata 将结果的第三类 (prog=='honors ') 视为 基本结果 ，并从其余参数中减去其所有参数。在 Python 中，您可以通过 运行
print(mul_lr.intercept_ - mul_lr.intercept_[-1])
print((mul_lr.coef_ - mul_lr.coef_[-1])[:, reorder])
这给了你
[9.53746174 5.0549659 0. ]
[[ 0.57899547 -0.4976741 -0.04032136 -0.0459539 -0.09877841]
[-0.25622441 -0.74825829 -0.02334893 -0.01148954 -0.05993533]
[ 0. 0. 0. 0. 0. ]]
您现在可以看到参数现在接近 Stata 给出的参数：
- 截取 Python 中的 (9.53, 5.05) 与 Stata 中的 (9.54, 5.04)
- first-outcome 系数 (0.57, -0.49, ...) 对比 (0.61, -0.51, ...)
- second-outcome 系数 (-0.25, -0.74, ...) 对比 (-0.33, -0.86, ...)
你能看出规律吗？在 sklearn 中，斜率系数比在 Stata 中更小（接近于零）。这不是意外！
3。处理正则化
发生这种情况是因为 sklearn 有意 将斜率系数缩小到 0，方法是在系数上对其最大化的似然函数添加二次惩罚。这使得估计有偏差但更稳定，即使在严重的多重共线性的情况下也是如此。用贝叶斯术语来说，这种正则化对应于所有系数的 zero-mean 高斯先验。您可以了解有关正则化的更多信息 in the wiki.
在sklearn中，这个二次惩罚由正C参数控制：它越小，你得到的正则化程度越高。您可以将其视为每个斜率系数的先验方差。默认值是 C=1，但是你可以把它变大，比如 C=1000000，这意味着几乎没有正则化。在这种情况下，输出几乎与 Stata:
mul_lr2 = linear_model.LogisticRegression(
multi_class='multinomial', solver='newton-cg', C=1000000
).fit(train_x, train_y)
print(mul_lr2.intercept_ - mul_lr2.intercept_[-1])
print((mul_lr2.coef_ - mul_lr2.coef_[-1])[:, reorder])
这给了你
[9.54412644 5.04126452 0. ]
[[ 0.61978951 -0.51320481 -0.04053013 -0.0459711 -0.09904948]
[-0.33508605 -0.86869799 -0.02262518 -0.01161839 -0.05913068]
[ 0. 0. 0. 0. 0. ]]
结果仍然略有不同（如小数点后 5 位），但正则化更少时，差异填充会进一步缩小。