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

Question

我正在尝试使用 bootstrap 计算逻辑回归系数的方差，并且我正在使用 scikit-learn 和 statsmodels 来比较结果。我正在使用默认数据集 ISLR 网站，可以在 zip forlder here 或 here 中以纯 csv 文件形式找到。 我正在使用以下代码来执行引导程序：

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

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 生成的值更接近真实值，而且对于不同的运行，scikit-learn 系数似乎差异很大从实际值。你能解释一下为什么会这样吗？我觉得这很令人困惑，因为我希望对于相同的数据集，结果应该是相同的（至少略有不同）。然而，在这种情况下，结果差异很大，这让我相信我运行 sk-learn 版本的方式有问题。非常感谢任何帮助，我非常乐意提供更多说明。

Answer 1

虽然您将 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

虽然现在的系数相似，但它表明您的数据集存在更大的问题，类似于question

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]])