sklearn과 statsmodels 비교¶

import pandas as pd

# MASS 패키지의 Boston 자료 (Rdatasets 미러). sklearn에서는 1.2판에 제거되었다.
url = ("https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/"
       "1dcc2bf5f955cc1224a3e1307256e1fe86b68dae/csv/MASS/Boston.csv")
boston = pd.read_csv(url)
df = boston.rename(columns={"rm": "RM", "lstat": "LSTAT",
                            "ptratio": "PTRATIO", "medv": "PRICE"})
df = df[["RM", "LSTAT", "PTRATIO", "PRICE"]]

print(df.describe().round(3).to_string())

            RM    LSTAT  PTRATIO    PRICE
count  506.000  506.000  506.000  506.000
mean     6.285   12.653   18.456   22.533
std      0.703    7.141    2.165    9.197
min      3.561    1.730   12.600    5.000
25%      5.885    6.950   17.400   17.025
50%      6.208   11.360   19.050   21.200
75%      6.624   16.955   20.200   25.000
max      8.780   37.970   22.000   50.000

from scipy import stats

print(f"PRICE 의 최댓값 = {df['PRICE'].max()},  정확히 50 인 건수 = {(df['PRICE'] == 50).sum()},  45 이상 = {(df['PRICE'] >= 45).sum()}")
for c in ["RM", "LSTAT", "PTRATIO", "PRICE"]:
    print(f"{c:8s} 평균 {df[c].mean():7.3f}  중앙값 {df[c].median():7.3f}  "
          f"평균-중앙값 {df[c].mean() - df[c].median():+7.3f}  왜도 {stats.skew(df[c]):+.3f}")
print(f"PTRATIO 의 서로 다른 값 = {df['PTRATIO'].nunique()}개,  최빈값 20.2 가 {(df['PTRATIO'] == 20.2).sum()}건")

PRICE 의 최댓값 = 50.0,  정확히 50 인 건수 = 16,  45 이상 = 22
RM       평균   6.285  중앙값   6.208  평균-중앙값  +0.076  왜도 +0.402
LSTAT    평균  12.653  중앙값  11.360  평균-중앙값  +1.293  왜도 +0.904
PTRATIO  평균  18.456  중앙값  19.050  평균-중앙값  -0.594  왜도 -0.800
PRICE    평균  22.533  중앙값  21.200  평균-중앙값  +1.333  왜도 +1.105
PTRATIO 의 서로 다른 값 = 46개,  최빈값 20.2 가 140건

import statsmodels.api as sm
import pandas as pd

# statsmodels 는 추론이 목적일 때 쓴다. 표준오차·p-값·신뢰구간이
# 출력표에 한꺼번에 나온다.
X = sm.add_constant(df[['RM', 'LSTAT', 'PTRATIO']])
y = df['PRICE']

def print_summary(res):
    """summary()에서 실행 날짜와 시각만 지우고 인쇄한다(재현 가능한 출력을 위해)."""
    lines = []
    for line in str(res.summary()).split("\n"):
        if line.startswith(("Date:", "Time:")):
            lines.append(line[:19].ljust(38) + line[38:])
        else:
            lines.append(line)
    print("\n".join(lines))


model = sm.OLS(y, X).fit()
print_summary(model)

# 필요한 값은 속성으로 바로 꺼낼 수 있다.
print(f"R²: {model.rsquared:.4f}")
print(f"Adj. R²: {model.rsquared_adj:.4f}")
print(f"AIC: {model.aic:.2f}")
print(f"BIC: {model.bic:.2f}")
print(model.pvalues)
print(model.conf_int())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                  PRICE   R-squared:                       0.679
Model:                            OLS   Adj. R-squared:                  0.677
Method:                 Least Squares   F-statistic:                     353.3
Date:                                   Prob (F-statistic):          2.69e-123
Time:                                   Log-Likelihood:                -1553.0
No. Observations:                 506   AIC:                             3114.
Df Residuals:                     502   BIC:                             3131.
Df Model:                           3                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         18.5671      3.913      4.745      0.000      10.879      26.255
RM             4.5154      0.426     10.603      0.000       3.679       5.352
LSTAT         -0.5718      0.042    -13.540      0.000      -0.655      -0.489
PTRATIO       -0.9307      0.118     -7.911      0.000      -1.162      -0.700
==============================================================================
Omnibus:                      202.072   Durbin-Watson:                   0.901
Prob(Omnibus):                  0.000   Jarque-Bera (JB):             1022.153
Skew:                           1.700   Prob(JB):                    1.10e-222
Kurtosis:                       9.076   Cond. No.                         402.
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
R²: 0.6786
Adj. R²: 0.6767
AIC: 3114.10
BIC: 3131.00
const      2.725808e-06
RM         7.734793e-24
LSTAT      7.944208e-36
PTRATIO    1.644660e-14
dtype: float64
                 0          1
const    10.878841  26.255382
RM        3.678711   5.352131
LSTAT    -0.654775  -0.488836
PTRATIO  -1.161877  -0.699568

import numpy as np
from statsmodels.stats.stattools import durbin_watson

e = np.asarray(model.resid)
dw_hand = (np.diff(e) ** 2).sum() / (e ** 2).sum()
rho = np.corrcoef(e[:-1], e[1:])[0, 1]
print(f"DW 손계산 = {dw_hand:.6f},  durbin_watson() = {durbin_watson(e):.6f}")
print(f"잔차의 1차 자기상관 rho = {rho:.6f},  2(1 - rho) = {2 * (1 - rho):.6f}")
print(f"분산 부풀림 (1+rho)/(1-rho) = {(1 + rho) / (1 - rho):.4f},  "
      f"표준오차 쪽은 그 제곱근 {np.sqrt((1 + rho) / (1 - rho)):.4f}")

# 자기상관에 강건한 표준오차와 견준다.
hac = sm.OLS(y, X).fit(cov_type='HAC', cov_kwds={'maxlags': 10})
tab = pd.DataFrame({'nonrobust': model.bse, 'HAC(10)': hac.bse})
tab['비'] = tab['HAC(10)'] / tab['nonrobust']
print(tab.round(4).to_string())
print("HAC 의 t 값:", hac.tvalues.round(3).tolist())
print("HAC 의 p 값:", ["%.2e" % v for v in hac.pvalues])

DW 손계산 = 0.901236,  durbin_watson() = 0.901236
잔차의 1차 자기상관 rho = 0.546931,  2(1 - rho) = 0.906137
분산 부풀림 (1+rho)/(1-rho) = 3.4143,  표준오차 쪽은 그 제곱근 1.8478
         nonrobust  HAC(10)       비
const       3.9132   9.9983  2.5550
RM          0.4259   1.4322  3.3629
LSTAT       0.0422   0.1108  2.6244
PTRATIO     0.1177   0.1682  1.4299
HAC 의 t 값: [1.857, 3.153, -5.159, -5.532]
HAC 의 p 값: ['6.33e-02', '1.62e-03', '2.48e-07', '3.16e-08']

from sklearn.linear_model import LinearRegression
from sklearn.model_selection import cross_val_score
from sklearn.metrics import r2_score, mean_squared_error
import numpy as np

# sklearn 은 예측이 목적일 때 쓴다. 표준오차나 p-값은 아예 제공하지 않는
# 대신, 교차검증·파이프라인·정규화가 한 틀로 묶여 있다.
X = df[['RM', 'LSTAT', 'PTRATIO']]
y = df['PRICE']

model = LinearRegression()
model.fit(X, y)

y_pred = model.predict(X)
print(f"R²: {r2_score(y, y_pred):.4f}")
print(f"RMSE: {np.sqrt(mean_squared_error(y, y_pred)):.4f}")
print(f"Coefficients: {model.coef_}")
print(f"Intercept: {model.intercept_:.4f}")

# 훈련자료에서 잰 R^2 는 언제나 후하다. 교차검증 값이 실제 성능에 가깝다.
cv_scores = cross_val_score(model, X, y, cv=5, scoring='r2')
print(f"CV R² (mean ± std): {cv_scores.mean():.4f} ± {cv_scores.std():.4f}")

R²: 0.6786
RMSE: 5.2087
Coefficients: [ 4.51542094 -0.57180569 -0.93072256]
Intercept: 18.5671
CV R² (mean ± std): 0.4300 ± 0.2997

from sklearn.model_selection import KFold

sm_fit = sm.OLS(df['PRICE'], sm.add_constant(df[['RM', 'LSTAT', 'PTRATIO']])).fit()
print(f"계수의 최대 차이   = {np.abs(np.r_[model.intercept_, model.coef_] - sm_fit.params.values).max():.3e}")
print(f"적합값의 최대 차이 = {np.abs(model.predict(X) - sm_fit.fittedvalues).max():.3e}")
print(f"훈련 R^2 = {r2_score(y, y_pred):.4f},  낙관의 어림값 2p/n = {2 * 4 / len(df):.4f}")

print("섞지 않은 CV:", cv_scores.round(4), f" 평균 {cv_scores.mean():.4f}  표준편차 {cv_scores.std():.4f}")
shuffled = cross_val_score(LinearRegression(), X, y,
                           cv=KFold(5, shuffle=True, random_state=0), scoring='r2')
print("섞은 CV     :", shuffled.round(4), f" 평균 {shuffled.mean():.4f}  표준편차 {shuffled.std():.4f}")
for i, (_, te) in enumerate(KFold(5).split(X)):
    print(f"  겹 {i}: 색인 {te.min():3d}~{te.max():3d}  PRICE 평균 {y.iloc[te].mean():.2f}")

계수의 최대 차이   = 7.105e-14
적합값의 최대 차이 = 8.660e-14
훈련 R^2 = 0.6786,  낙관의 어림값 2p/n = 0.0158
섞지 않은 CV: [ 0.7269  0.7149  0.5482  0.1819 -0.0219]  평균 0.4300  표준편차 0.2997
섞은 CV     : [0.4882 0.7613 0.6543 0.6005 0.7952]  평균 0.6599  표준편차 0.1111
  겹 0: 색인   0~101  PRICE 평균 22.40
  겹 1: 색인 102~202  PRICE 평균 24.42
  겹 2: 색인 203~303  PRICE 평균 29.70
  겹 3: 색인 304~404  PRICE 평균 20.04
  겹 4: 색인 405~505  PRICE 평균 16.10

# 실무에서는 둘을 함께 쓴다. 먼저 statsmodels 로 무엇이 유의한지,
# 가정이 지켜지는지를 살핀다.
import pandas as pd
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor

X_sm = sm.add_constant(df[['RM', 'LSTAT', 'PTRATIO']])
model_sm = sm.OLS(df['PRICE'], X_sm).fit()
print_summary(model_sm)

# 다중공선성 점검. 상수항은 셈에서 뺀다.
vif = pd.DataFrame({
    'Feature': X_sm.columns[1:],
    'VIF': [variance_inflation_factor(X_sm.values, i) for i in range(1, X_sm.shape[1])]
})
print(vif)

# 그다음 sklearn 으로 예측 파이프라인을 세운다. 표준화와 능형회귀를
# 하나로 묶으면, 교차검증의 각 겹에서 표준화가 훈련 부분만 보고 이뤄진다.
# 이렇게 해야 검증자료의 정보가 새어 들어가지 않는다.
from sklearn.linear_model import Ridge
from sklearn.model_selection import cross_val_score
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline

pipe = Pipeline([
    ('scaler', StandardScaler()),
    ('ridge', Ridge(alpha=1.0))
])

X_sk = df[['RM', 'LSTAT', 'PTRATIO']]
cv_scores = cross_val_score(pipe, X_sk, df['PRICE'], cv=5, scoring='r2')
print(f"Ridge CV R²: {cv_scores.mean():.4f} ± {cv_scores.std():.4f}")

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                  PRICE   R-squared:                       0.679
Model:                            OLS   Adj. R-squared:                  0.677
Method:                 Least Squares   F-statistic:                     353.3
Date:                                   Prob (F-statistic):          2.69e-123
Time:                                   Log-Likelihood:                -1553.0
No. Observations:                 506   AIC:                             3114.
Df Residuals:                     502   BIC:                             3131.
Df Model:                           3                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         18.5671      3.913      4.745      0.000      10.879      26.255
RM             4.5154      0.426     10.603      0.000       3.679       5.352
LSTAT         -0.5718      0.042    -13.540      0.000      -0.655      -0.489
PTRATIO       -0.9307      0.118     -7.911      0.000      -1.162      -0.700
==============================================================================
Omnibus:                      202.072   Durbin-Watson:                   0.901
Prob(Omnibus):                  0.000   Jarque-Bera (JB):             1022.153
Skew:                           1.700   Prob(JB):                    1.10e-222
Kurtosis:                       9.076   Cond. No.                         402.
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
   Feature       VIF
0       RM  1.653419
1    LSTAT  1.679425
2  PTRATIO  1.198101
Ridge CV R²: 0.4304 ± 0.2994

import numpy as np

Xv = X_sm.values
for i in range(1, X_sm.shape[1]):
    others = [j for j in range(X_sm.shape[1]) if j != i]
    aux = sm.OLS(Xv[:, i], Xv[:, others]).fit()
    print(f"  {X_sm.columns[i]:8s} VIF = {variance_inflation_factor(Xv, i):.6f},  "
          f"1/(1 - R^2) = {1 / (1 - aux.rsquared):.6f}")

print(f"조건수 그대로         = {model_sm.condition_number:.2f}")
print(f"열별 노름             = {np.linalg.norm(Xv, axis=0).round(2)}")
print(f"열 노름을 1 로 맞추면  = {np.linalg.cond(Xv / np.linalg.norm(Xv, axis=0)):.3f}")
Z = (Xv[:, 1:] - Xv[:, 1:].mean(0)) / Xv[:, 1:].std(0)
print(f"중심화·표준화 뒤       = {np.linalg.cond(np.column_stack([np.ones(len(Z)), Z])):.3f}")

d = np.linalg.svd(Z[:405], compute_uv=False)
print(f"표준화 열의 특이값(훈련 405행) = {d.round(1)},  d^2 = {(d ** 2).round(0)}")
print(f"릿지의 축소율 d^2/(d^2 + 1)    = {(d ** 2 / (d ** 2 + 1)).round(5)}")
ols_cv = cross_val_score(LinearRegression(), X_sk, df['PRICE'], cv=5, scoring='r2')
print(f"Ridge CV {cv_scores.mean():.4f} +- {cv_scores.std():.4f},  "
      f"OLS CV {ols_cv.mean():.4f} +- {ols_cv.std():.4f}")

  RM       VIF = 1.653419,  1/(1 - R^2) = 1.653419
  LSTAT    VIF = 1.679425,  1/(1 - R^2) = 1.679425
  PTRATIO  VIF = 1.198101,  1/(1 - R^2) = 1.198101
조건수 그대로         = 401.61
열별 노름             = [ 22.49 142.25 326.75 417.99]
열 노름을 1 로 맞추면  = 39.888
중심화·표준화 뒤       = 2.224
표준화 열의 특이값(훈련 405행) = [28.3 17.9 12.1],  d^2 = [803. 320. 146.]
릿지의 축소율 d^2/(d^2 + 1)    = [0.99876 0.99688 0.99319]
Ridge CV 0.4304 +- 0.2994,  OLS CV 0.4300 +- 0.2997