来自 sk_learn 的 LogisticRegression 和来自 statsmodels.formula.api return 的 smf.logit() 不同的结果

LogisticRegression from sk_learn and smf.logit() from statsmodels.formula.api return different results

我正在尝试使用 bootstrap 计算逻辑回归系数的方差,并且我正在使用 scikit-learn 和 statsmodels 来比较结果。我正在使用来自 ISLR website which can be found in the zip forlder here or here 作为普通的 csv 文件。 我正在使用以下代码执行 bootstrap:

导入数据集并创建响应变量

default_df = pd.read_csv("./Machine-Learning-Books-With-Python/Introduction to Statistical Learning/data/default.csv")

default_df['default_01'] = np.where(default_df.default == 'Yes', 1, 0)

接下来,我正在定义 boot 函数,它将处理我的数据集的随机抽样:

def boot(X, bootSample_size=None):
    '''
    Sampling observations from a dataframe

    Parameters
    ------------
    X : pandas dataframe
        Data to be resampled

    bootSample_size: int, optional
        Dimension of the bootstrapped samples

    Returns
    ------------
    bootSample_X : pandas dataframe
                   Resampled data

    Examples
    ----------
    To resample data from the X dataframe:
        >> boot(X)
    The resampled data will have length equal to len(X).

    To resample data from the X dataframe in order to have length 5:
        >> boot(X,5)

    References
    ------------
    http://nbviewer.jupyter.org/gist/aflaxman/6871948

    '''
    #assign default size if non-specified
    if bootSample_size == None:
        bootSample_size = len(X)

    #create random integers to use as indices for bootstrap sample based on original data
    bootSample_i = (np.random.rand(bootSample_size)*len(X)).astype(int)
    bootSample_i = np.array(bootSample_i)
    bootSample_X = X.iloc[bootSample_i]

    return bootSample_X

最后,我定义了两个将执行逻辑回归并提取参数的函数:

使用 statsmodels

def boot_fn(data):
    lr = smf.logit(formula='default_01 ~ income + balance', data=data).fit(disp=0)
    return lr.params

使用 scikit-learn

def boot_fn2(data):
    X = data[['income', 'balance']]
    y = data.default_01
    logit = LogisticRegression(C = 1e9)
    logit.fit(X, y)
    return logit.coef_

最后循环 运行 函数 100 次并存储结果:

coef_sk = []
coef_sm = []
for _ in np.arange(100):
    data = boot(default_df)
    coef_sk.append(boot_fn2(data))
    coef_sm.append(boot_fn(data))

取 coef_sk (scikit-learn) 和 coef_sm (statsmodels) 的平均值,我发现使用 statsmodels 生成的值更接近实际值,并且对于不同的 运行s scikit-learn 系数似乎与实际值有很大差异。你能解释一下为什么会这样吗?我觉得这很混乱,因为我希望对于相同的数据集,结果应该是相同的(至少略有不同)。然而,在这种情况下,结果相差很大,这让我相信我 运行ning sk-learn 版本的方式有问题。非常感谢任何帮助,我非常乐意提供更多说明。

虽然您将 C 参数设置为高以最小化,但 sklearn 默认使用 lbfgs 求解器来找到您的最佳参数,而 statsmodels 使用 newton .

您可以尝试这样做以获得相似的系数:

def boot_fn2(data):
    X = data[['income', 'balance']]
    y = data.default_01
    logit = LogisticRegression(penalty="none",max_iter=1000,solver = "newton-cg")
    logit.fit(X, y)
    return logit.coef_

如果我运行这个有上面的功能:

coef_sk = []
coef_sm = []
for _ in np.arange(50):
    data = boot(default_df)
    coef_sk.append(boot_fn2(data))
    coef_sm.append(boot_fn(data))

你会立即看到它抛出了很多关于无法收敛的警告:

LineSearchWarning: The line search algorithm did not converge

尽管现在系数相似,但它表明您的数据集存在更大的问题,类似于此

np.array(coef_sm)[:,1:].mean(axis=0)
array([2.14570133e-05, 5.68280785e-03])

np.array(coef_sk).mean(axis=0)
array([[2.14352318e-05, 5.68116402e-03]])

你的因变量非常大,这给 sklearn 中可用的优化方法带来了问题。如果你想解释系数,你可以缩小你的两个因变量:

default_df[['balance','income']] = default_df[['balance','income']]/100

否则,首先缩放自变量并应用回归始终是一个好习惯:

from sklearn.preprocessing import StandardScaler
default_df[['balance','income']] = StandardScaler().fit_transform(default_df[['balance','income']])

def boot_fn(data):
    lr = smf.logit(formula='default_01 ~ income + balance', data=data).fit(disp=0)
    return lr.params
def boot_fn2(data):
    X = data[['income', 'balance']]
    y = data.default_01
    logit = LogisticRegression(penalty="none")
    logit.fit(X, y)
    return logit.coef_

coef_sk = []
coef_sm = []
for _ in np.arange(50):
    data = boot(default_df)
    #print(data.default_01.mean())
    coef_sk.append(boot_fn2(data))
    coef_sm.append(boot_fn(data))

现在您会看到系数相似:

np.array(coef_sm)[:,1:].mean(axis=0)    
array([0.26517582, 2.71598194])

np.array(coef_sk).mean(axis=0)
array([[0.26517504, 2.71598548]])