statsmodels 최소제곱 인터페이스¶

import numpy as np
import statsmodels.api as sm

# 참 계수가 [3.0, 1.5], 절편이 2.0 인 자료다.
np.random.seed(42)
n = 100
X = np.random.randn(n, 2)
beta_true = np.array([3.0, 1.5])
y = X @ beta_true + 2.0 + np.random.randn(n) * 0.5

# statsmodels 는 절편을 자동으로 넣지 않는다. 이 줄을 빠뜨리면 원점을
# 지나는 회귀가 되므로 sklearn 과의 가장 흔한 차이가 여기서 생긴다.
X_with_const = sm.add_constant(X)

# OLS(y, X) 순서다. sklearn 의 fit(X, y) 와 반대이니 헷갈리기 쉽다.
model = sm.OLS(y, X_with_const)
results = model.fit()

def print_summary(res):
    """summary()에서 실행 날짜와 시각만 지우고 인쇄한다.

    statsmodels의 summary()는 표 머리에 Date와 Time을 함께 찍는다.
    그대로 두면 실행할 때마다 출력이 달라져 문서에 싣기 어렵다.
    같은 줄에 있는 다른 값(Prob (F-statistic), Log-Likelihood)은 남긴다.
    """
    lines = []
    for line in str(res.summary()).split("\n"):
        if line.startswith(("Date:", "Time:")):
            label, rest = line[:19], line[19:]
            lines.append(label.ljust(19) + " " * 19 + rest[19:])
        else:
            lines.append(line)
    print("\n".join(lines))

print_summary(results)

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.971
Model:                            OLS   Adj. R-squared:                  0.970
Method:                 Least Squares   F-statistic:                     1603.
Date:                                   Prob (F-statistic):           4.97e-75
Time:                                   Log-Likelihood:                -77.799
No. Observations:                 100   AIC:                             161.6
Df Residuals:                      97   BIC:                             169.4
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          2.0464      0.054     37.884      0.000       1.939       2.154
x1             3.0954      0.063     49.281      0.000       2.971       3.220
x2             1.4139      0.054     26.258      0.000       1.307       1.521
==============================================================================
Omnibus:                        4.136   Durbin-Watson:                   2.212
Prob(Omnibus):                  0.126   Jarque-Bera (JB):                3.956
Skew:                           0.266   Prob(JB):                        0.138
Kurtosis:                       3.817   Cond. No.                         1.23
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

from scipy import stats

p = X_with_const.shape[1]
e = results.resid
rss = (e ** 2).sum()
s = np.sqrt(rss / (n - p))
XtX_inv = np.linalg.inv(X_with_const.T @ X_with_const)

se = s * np.sqrt(np.diag(XtX_inv))
print(f"s = sqrt(RSS/(n-p)) = {s:.6f}")
print("SE  손계산:", se.round(6), "  statsmodels:", results.bse.round(6))
t = results.params / se
print("t   손계산:", t.round(4), "  statsmodels:", results.tvalues.round(4))
t_crit = stats.t.ppf(0.975, n - p)
ci = np.column_stack([results.params - t_crit * se, results.params + t_crit * se])
print(f"t 임계값 t_(0.975, {n - p}) = {t_crit:.5f},  신뢰구간 최대 차이 = {np.abs(ci - results.conf_int()).max():.3e}")

R2 = results.rsquared
print(f"F   손계산 = {(R2 / (p - 1)) / ((1 - R2) / (n - p)):.4f}   statsmodels = {results.fvalue:.4f}")
print(f"adj 손계산 = {1 - (1 - R2) * (n - 1) / (n - p):.10f}   statsmodels = {results.rsquared_adj:.10f}")
print(f"logL 손계산 = {-n / 2 * (np.log(2 * np.pi) + np.log(rss / n) + 1):.6f}   statsmodels = {results.llf:.6f}")
print(f"AIC 손계산 = {-2 * results.llf + 2 * p:.6f}   statsmodels = {results.aic:.6f}")
print(f"BIC 손계산 = {-2 * results.llf + p * np.log(n):.6f}   statsmodels = {results.bic:.6f}")
print(f"조건수 손계산 = {np.linalg.cond(X_with_const):.4f}   statsmodels = {results.condition_number:.4f}")

s = sqrt(RSS/(n-p)) = 0.534875
SE  손계산: [0.054017 0.06281  0.053847]   statsmodels: [0.054017 0.06281  0.053847]
t   손계산: [37.8844 49.2813 26.2585]   statsmodels: [37.8844 49.2813 26.2585]
t 임계값 t_(0.975, 97) = 1.98472,  신뢰구간 최대 차이 = 2.220e-16
F   손계산 = 1602.5775   statsmodels = 1602.5775
adj 손계산 = 0.9700195786   statsmodels = 0.9700195786
logL 손계산 = -77.798698   statsmodels = -77.798698
AIC 손계산 = 161.597396   statsmodels = 161.597396
BIC 손계산 = 169.412907   statsmodels = 169.412907
조건수 손계산 = 1.2339   statsmodels = 1.2339

import pandas as pd
import statsmodels.formula.api as smf

# 수식 API 는 R 의 문법을 따른다. 이쪽에서는 절편이 자동으로 들어가고,
# 범주형 변수도 알아서 가변수로 바뀐다.
df = pd.DataFrame({
    'y': y,
    'x1': X[:, 0],
    'x2': X[:, 1]
})

results_formula = smf.ols('y ~ x1 + x2', data=df).fit()
print_summary(results_formula)

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.971
Model:                            OLS   Adj. R-squared:                  0.970
Method:                 Least Squares   F-statistic:                     1603.
Date:                                   Prob (F-statistic):           4.97e-75
Time:                                   Log-Likelihood:                -77.799
No. Observations:                 100   AIC:                             161.6
Df Residuals:                      97   BIC:                             169.4
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept      2.0464      0.054     37.884      0.000       1.939       2.154
x1             3.0954      0.063     49.281      0.000       2.971       3.220
x2             1.4139      0.054     26.258      0.000       1.307       1.521
==============================================================================
Omnibus:                        4.136   Durbin-Watson:                   2.212
Prob(Omnibus):                  0.126   Jarque-Bera (JB):                3.956
Skew:                           0.266   Prob(JB):                        0.138
Kurtosis:                       3.817   Cond. No.                         1.23
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

print(f"설계행렬의 최대 차이 = {np.abs(np.asarray(results_formula.model.exog) - X_with_const).max():.3e}")
print("열 이름:", results_formula.model.exog_names, "대", results.model.exog_names)
for attr in ['params', 'bse', 'tvalues', 'pvalues']:
    d = np.abs(np.asarray(getattr(results_formula, attr)) - np.asarray(getattr(results, attr))).max()
    print(f"  {attr:9s} 최대 차이 = {d:.3e}")
for attr in ['rsquared', 'fvalue', 'llf', 'aic', 'bic']:
    d = abs(getattr(results_formula, attr) - getattr(results, attr))
    print(f"  {attr:9s} 차이 = {d:.3e}")

# 절편을 지우면 R^2 의 정의가 바뀐다.
no_int = smf.ols('y ~ x1 + x2 - 1', data=df).fit()
print(f"'- 1' 모형: 설명변수 {no_int.model.exog_names},  k_constant = {no_int.k_constant}")
print(f"  statsmodels 가 찍는 R^2        = {no_int.rsquared:.6f}")
print(f"  1 - sum e^2 / sum y^2          = {1 - (no_int.resid ** 2).sum() / (y ** 2).sum():.6f}")
print(f"  같은 모형을 중심화해 다시 계산  = {1 - (no_int.resid ** 2).sum() / ((y - y.mean()) ** 2).sum():.6f}")
print(f"  잔차의 합: 절편 있음 {results.resid.sum():.3e},  절편 없음 {no_int.resid.sum():.4f}")

설계행렬의 최대 차이 = 0.000e+00
열 이름: ['Intercept', 'x1', 'x2'] 대 ['const', 'x1', 'x2']
  params    최대 차이 = 0.000e+00
  bse       최대 차이 = 0.000e+00
  tvalues   최대 차이 = 0.000e+00
  pvalues   최대 차이 = 0.000e+00
  rsquared  차이 = 0.000e+00
  fvalue    차이 = 0.000e+00
  llf       차이 = 0.000e+00
  aic       차이 = 0.000e+00
  bic       차이 = 0.000e+00
'- 1' 모형: 설명변수 ['x1', 'x2'],  k_constant = 0
  statsmodels 가 찍는 R^2        = 0.648290
  1 - sum e^2 / sum y^2          = 0.648290
  같은 모형을 중심화해 다시 계산  = 0.535991
  잔차의 합: 절편 있음 2.204e-14,  절편 없음 200.6488

# 적합 결과에서 꺼낼 수 있는 것들을 한자리에 모았다.
print("Coefficients:", results.params)

# 표준오차
print("Standard errors:", results.bse)

# p-값
print("P-values:", results.pvalues)

# 신뢰구간
print("95% CI:\n", results.conf_int(alpha=0.05))

# 결정계수와 수정결정계수
print("R-squared:", results.rsquared)
print("Adjusted R-squared:", results.rsquared_adj)

# 잔차
residuals = results.resid

# 적합값
fitted = results.fittedvalues

# AIC and BIC
print("AIC:", results.aic)
print("BIC:", results.bic)

Coefficients: [2.04639669 3.09536017 1.41392895]
Standard errors: [0.05401682 0.06281005 0.05384654]
P-values: [6.18804800e-60 1.81607514e-70 7.11568582e-46]
95% CI:
 [[1.93918825 2.15360513]
 [2.9706996  3.22002073]
 [1.30705847 1.52079942]]
R-squared: 0.9706252436321817
Adjusted R-squared: 0.9700195785524329
AIC: 161.5973963941535
BIC: 169.41290695211777

print(f"신뢰구간 = params +- t_crit * bse  최대 차이 = "
      f"{np.abs(np.column_stack([results.params - t_crit * results.bse, results.params + t_crit * results.bse]) - results.conf_int()).max():.3e}")
pv = 2 * stats.t.sf(np.abs(results.tvalues), n - p)
print(f"p-값 = 2 * P(T_{n - p} > |t|)      최대 상대차이 = {np.abs(pv / results.pvalues - 1).max():.3e}")
print(f"adj R^2 = 1 - (1-R^2)(n-1)/(n-p) = {1 - (1 - results.rsquared) * (n - 1) / (n - p):.10f}")
print(f"잔차의 합 = {results.resid.sum():.3e},  적합값 + 잔차 = y 인가: "
      f"{np.abs(results.fittedvalues + results.resid - y).max():.3e}")

신뢰구간 = params +- t_crit * bse  최대 차이 = 0.000e+00
p-값 = 2 * P(T_97 > |t|)      최대 상대차이 = 0.000e+00
adj R^2 = 1 - (1-R^2)(n-1)/(n-p) = 0.9700195786
잔차의 합 = 2.204e-14,  적합값 + 잔차 = y 인가: 2.776e-17

# 진단 도구가 갖춰져 있다는 점이 statsmodels 를 쓰는 큰 이유다.
# sklearn 에는 이런 것이 아예 없다.
influence = results.get_influence()
cooks_d = influence.cooks_distance[0]
leverage = influence.hat_matrix_diag

# 등분산성 검정
from statsmodels.stats.diagnostic import het_breuschpagan
bp_stat, bp_pval, _, _ = het_breuschpagan(results.resid, results.model.exog)
print(f"Breusch-Pagan p-value: {bp_pval:.4f}")

# 잔차의 정규성 검정
from statsmodels.stats.stattools import jarque_bera
jb_stat, jb_pval, skew, kurtosis = jarque_bera(results.resid)
print(f"Jarque-Bera p-value: {jb_pval:.4f}")

# 다중공선성 점검. 상수항 열은 건너뛴다.
from statsmodels.stats.outliers_influence import variance_inflation_factor
for i in range(1, X_with_const.shape[1]):
    vif = variance_inflation_factor(X_with_const, i)
    print(f"VIF for variable {i}: {vif:.2f}")

Breusch-Pagan p-value: 0.7942
Jarque-Bera p-value: 0.1383
VIF for variable 1: 1.00
VIF for variable 2: 1.00

bp_stat, bp_pval, _, _ = het_breuschpagan(results.resid, X_with_const)
aux = sm.OLS(results.resid ** 2, X_with_const).fit()
print(f"BP 통계량 = {bp_stat:.6f},  n * R^2(보조회귀) = {n * aux.rsquared:.6f}")
print(f"BP p-값   = {bp_pval:.6f},  chi2(2) 의 꼬리 = {stats.chi2.sf(n * aux.rsquared, 2):.6f}")

jb_stat, jb_pval, skew, kurtosis = jarque_bera(results.resid)
print(f"왜도 = {skew:.6f},  첨도 = {kurtosis:.6f}")
print(f"JB 통계량 = {jb_stat:.6f},  손계산 = {n / 6 * (skew ** 2 + (kurtosis - 3) ** 2 / 4):.6f}")

r12 = np.corrcoef(X[:, 0], X[:, 1])[0, 1]
for i in [1, 2]:
    others = [j for j in range(p) if j != i]
    aux2 = sm.OLS(X_with_const[:, i], X_with_const[:, others]).fit()
    print(f"VIF{i} = {variance_inflation_factor(X_with_const, i):.6f},  "
          f"1/(1 - R^2) = {1 / (1 - aux2.rsquared):.6f},  1/(1 - r^2) = {1 / (1 - r12 ** 2):.6f}")

inf = results.get_influence()
h = inf.hat_matrix_diag
cooks_d2 = inf.cooks_distance[0]
s2 = rss / (n - p)
cooks_hand = results.resid ** 2 * h / (p * s2 * (1 - h) ** 2)
print(f"Cook 거리 최대 차이 = {np.abs(cooks_d2 - cooks_hand).max():.3e}")
print(f"레버리지의 합 = {h.sum():.6f}  (= p = {p}),  최대 Cook D = {cooks_d2.max():.4f} (4/n = {4 / n:.3f})")

BP 통계량 = 0.460940,  n * R^2(보조회귀) = 0.460940
BP p-값   = 0.794160,  chi2(2) 의 꼬리 = 0.794160
왜도 = 0.265727,  첨도 = 3.816738
JB 통계량 = 3.956269,  손계산 = 3.956269
VIF1 = 1.001046,  1/(1 - R^2) = 1.001046,  1/(1 - r^2) = 1.001046
VIF2 = 1.001046,  1/(1 - R^2) = 1.001046,  1/(1 - r^2) = 1.001046
Cook 거리 최대 차이 = 5.551e-17
레버리지의 합 = 3.000000  (= p = 3),  최대 Cook D = 0.1961 (4/n = 0.040)