无法在Sklearn中复现Statsmodels逻辑回归的AIC等指标，求排查

我在Sklearn中手动计算逻辑回归的AIC、R²、调整后R²、对数似然等指标时，结果和Statsmodels的输出完全不匹配，请帮忙排查代码错误。

原始代码 #

# Import
import statsmodels.api as sm
from statsmodels.tools import add_constant

import numpy as np
import pandas as pd

from math import log

from sklearn import linear_model
from sklearn.linear_model import LogisticRegressionCV
from sklearn.metrics import log_loss 
from sklearn.metrics import mean_squared_error

# Dataset
n = 100
data = {'CO2': np.random.normal(14, 3, n),
        'NO2': np.random.normal(96, 2, n),
        'Outcome': np.random.choice([0,1], 
                                  size=100, 
                                  p=[0.7, 0.3])}
df = pd.DataFrame(data=data)
print(df.head())

# Split
X = df.drop(columns=['Outcome'])
Y = df[['Outcome']]
Y = Y.values.ravel()

# sklearn 
logreg = LogisticRegressionCV(penalty="l1", fit_intercept = True, 
                              solver="saga", 
                              cv=5, random_state=0, max_iter=1000, scoring='accuracy')
logreg.fit(X, Y)

# negative log-likelihood of a logistic model that returns y_pred probabilities for its training data y_true

yhat = logreg.predict(X)
Y_yhat_df = pd.DataFrame({'Y': Y, 'yhat': yhat})
# print(Y_yhat_df.head())

logloss = round(log_loss(Y, yhat), 4)
print('log loss: %s' % logloss)

loglike = -1 * logloss
print('log likelihood: %s' % loglike)

mse = mean_squared_error(Y, yhat)
print('MSE: %s' % round(mse, 2))

# Model parameters
k = X.shape[1]
print('k: %s' % k)

n = len(X)
print('n: %s' % n)

num_params = logreg.coef_.shape[-1]
print('Parameters: %s' % num_params)
# print(logreg.coef_)

print('R2: %s' % round(logreg.score(X, Y), 2))

# Adjusted R^2
R2adj = 1 - (1-logreg.score(X, Y))*(len(Y)-1)/(len(Y)-X.shape[1]-1)
print('R2 adj: %s' % round(R2adj, 2))

# 2k - 2(log(log likelihood))
# np.log is the natural log
aic = 2*num_params - 2 * (-1 * np.log(logloss))
print('AIC: %s' % round(aic,2))

bic = n * log(mse) + num_params * log(n)
print('BIC: %s' % round(bic, 2))

# Stats models

# Manually add intercept 
Xi = sm.add_constant(X)

sm_logreg = sm.Logit(Y, Xi)

res = sm_logreg.fit_regularized(method='l1')

print(res.aic)
print(res.summary())

原始输出对比 #

Sklearn输出 #

• log loss: 12.6153
• log likelihood: -12.6153
• MSE: 0.35
• k: 2
• n: 100
• Parameters: 2
• R2: 0.65
• R2 adj: 0.64
• AIC: 9.07 BIC: -95.77

Statsmodels输出 #

• AIC: 134.88348802189196
• Pseudo R-squ.: 0.004679（调整后R²）
• Log-Likelihood: -64.442

核心错误点分析 #

1. 对数似然与Log Loss计算错误 #

• Sklearn的log_loss默认返回平均负对数似然，而Statsmodels的对数似然是总负对数似然的相反数，直接取反平均log loss完全错误。
• 你传入了predict(X)得到的类别标签（0/1），但log_loss需要模型输出的概率值，应该用predict_proba(X)[:, 1]。

2. R²与调整R²的误用 #

• Sklearn的logreg.score(X,Y)返回的是分类准确率，不是线性回归的R²。逻辑回归是分类模型，必须用伪R²（如McFadden R²），不能套用线性回归的R²公式。

3. AIC公式错误 #

• AIC正确公式为AIC = 2*k - 2*loglike，其中k是总参数数（含截距），你漏掉了截距项，还错误地对log loss取自然对数，完全偏离定义。

4. BIC公式错误 #

• 分类模型的BIC公式为BIC = -2*loglike + k*log(n)，你用了线性回归的MSE相关公式，不适用逻辑回归。

5. 模型正则化参数不统一 #

• Sklearn和Statsmodels的L1正则化默认强度不同，导致模型参数不一致，指标自然无法匹配。

修正后的完整代码 #

# Import
import statsmodels.api as sm
from statsmodels.tools import add_constant
from statsmodels.discrete.discrete_model import Logit

import numpy as np
import pandas as pd

from sklearn.linear_model import LogisticRegressionCV
from sklearn.metrics import log_loss 

# Dataset
n = 100
np.random.seed(42)  # 设置随机种子保证可复现
data = {'CO2': np.random.normal(14, 3, n),
        'NO2': np.random.normal(96, 2, n),
        'Outcome': np.random.choice([0,1], 
                                  size=100, 
                                  p=[0.7, 0.3])}
df = pd.DataFrame(data=data)

# Split
X = df.drop(columns=['Outcome'])
Y = df['Outcome'].values.ravel()

# sklearn 
logreg = LogisticRegressionCV(penalty="l1", fit_intercept = True, 
                              solver="saga", 
                              cv=5, random_state=0, max_iter=1000, scoring='accuracy')
logreg.fit(X, Y)

# 用概率值计算总负对数似然
y_proba = logreg.predict_proba(X)[:, 1]
logloss_total = log_loss(Y, y_proba, normalize=False)
loglike = -logloss_total
print(f'对数似然: {round(loglike, 3)}')

# 计算McFadden伪R²
X_null = sm.add_constant(np.ones(n))
null_model = Logit(Y, X_null).fit(disp=False)
loglike_null = null_model.llf
mcfadden_r2 = 1 - (loglike / loglike_null)
print(f'McFadden伪R²: {round(mcfadden_r2, 4)}')

# 计算AIC和BIC（包含截距项）
k = logreg.coef_.shape[-1] + 1  
aic = 2*k - 2*loglike
print(f'AIC: {round(aic, 2)}')
bic = -2*loglike + k*np.log(n)
print(f'BIC: {round(bic, 2)}')

# Stats models
Xi = sm.add_constant(X)
sm_logreg = sm.Logit(Y, Xi)
res = sm_logreg.fit_regularized(method='l1', disp=False)

print('\nStatsmodels输出:')
print(f'AIC: {round(res.aic, 2)}')
print(f'对数似然: {round(res.llf, 3)}')
print(f'McFadden伪R²: {round(res.prsquared, 4)}')

内容的提问来源于stack exchange，提问作者Sandra T