多项式Logit模型：Python和Stata的结果不同

Question

多项式Logit模型：Python和Stata的结果不同

7

我正在尝试使用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-learn模块复制相同的结果。以下是代码：

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 是一个指示变量？

编辑：

数据文件链接

- Furqan Hashim

这是反向的。除非Python正确地指示了SES，否则会得到不同的结果。您的第二个问题仍然存在。 - Nick Cox

如果您想模仿Stata输出，您需要为ses==1和ses==2构建虚拟代码。请参见例如：https://pandas.pydata.org/pandas-docs/stable/generated/pandas.get_dummies.html - Richard Border

使用sklearn.preprocessing.LabelEncoder得到了相同的结果。尝试使用pandas get_dummies，但它会为三个类别创建三个单独的列，而不是创建一个编码为ses==1 ses==2的单个列。仍在努力复制相同的结果。 - Furqan Hashim

您附加的数据编码不同：ses具有“低”，“中”和“高”级别。请以相同的格式提供数据和代码！ - David Dale

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- David Dale · Accepted Answer

有几个问题会导致 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）
第一个结果系数为（0.57，-0.49，...），而Stata中为（0.61，-0.51，...）
第二个结果系数为（-0.25，-0.74，...），而Stata中为（-0.33，-0.86，...）

你能看到规律吗？在sklearn中，斜率系数比Stata小（接近于零）。这不是偶然的！

3. 处理正则化

这是因为sklearn有意通过在其最大化的似然函数中添加二次惩罚项来将斜率系数向0缩小。这使得估计结果具有偏差但更加稳定，即使在存在严重多重共线性的情况下也是如此。从贝叶斯的角度来看，这种正则化对应于所有系数的零均值高斯先验分布。您可以在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位小数处），但是使用更少的正则化，差异会进一步缩小。