ADP 실기 — 회귀분석, 적합보다 진단이 배점이다

LinearRegression().fit(X, y) 는 한 줄이다. 시험에서 그 한 줄에는 점수가 거의 없다. 회귀계수를 검정하고, 모형을 검정하고, 가정이 지켜졌는지 확인하는 것에 점수가 있다.

그래서 회귀 문제는 sklearn 이 아니라 statsmodels 로 푼다. sklearn 의 LinearRegression 은 p 값을 주지 않는다.

다중회귀 — 기본형

import statsmodels.api as sm
from statsmodels.formula.api import ols

model = ols('value ~ item + layer + item:layer', data=DF).fit()
print(model.summary())

print(sm.stats.anova_lm(model))          # F 검정 — 모형 전체의 유의성
print(model.t_test([0, 1, -1, 0]))       # T 검정 — 계수의 선형 결합

ols 의 문자열 수식이 편한 이유는 범주형을 알아서 더미로 만들어 주기 때문이다. C(gender) 로 감싸면 숫자로 들어 있는 열도 범주로 취급한다.

sklearn 으로 풀어야 한다면 더미화를 직접 한다.

df = pd.get_dummies(df, columns=['gender', 'city'], drop_first=True)

24회의 “광고비를 가변수화해서 다중선형회귀방정식을 만들고 회귀계수를 검정하기”가 정확히 이 두 줄이다. 가변수화 → 적합 → summary() 의 P>|t| 열을 읽는다.

다중공선성 — VIF

from statsmodels.stats.outliers_influence import variance_inflation_factor

X = df.drop('income', axis=1)
vif = pd.DataFrame()
vif['VIF Factor'] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
vif['features'] = X.columns
print(vif)

VIF 가 10 을 넘으면 다중공선성을 의심한다. 대처는 셋 중 하나다 — 변수를 빼거나, PCA 로 차원을 줄이거나, 릿지·라쏘로 벌점을 준다.

$\text{VIF}_i = \frac{1}{1 - R_i^2}$

$R_i^2$ 은 $i$ 번째 변수를 나머지 변수로 회귀했을 때의 결정계수다. 다른 변수들로 그 변수를 거의 설명할 수 있으면 VIF 가 치솟는다는 뜻이라, 정의만 봐도 왜 10 이 기준인지 짐작이 간다 ( $R^2 = 0.9$ 일 때 VIF 가 10 이다).

가정 다섯 가지를 검증한다

회귀분석 문제에서 가장 확실하게 점수가 나오는 부분이다. 순서대로 그리고 찍는다.

1. 선형성 — 잔차 대 적합값

y_pred = model.predict()
residuals = model.resid

sns.regplot(x=y_pred, y=residuals, lowess=True, line_kws={'color': 'red'})
plt.title('Residuals vs Fitted')
plt.show()

빨간 선이 0 근처에서 평평하면 선형성을 만족한다. 곡선이 보이면 변수를 변환하거나 다항항을 넣어야 한다.

2. 등분산성

from statsmodels.stats.diagnostic import het_breuschpagan

lm, lm_p, f, f_p = het_breuschpagan(model.resid, model.model.exog)
print(lm_p)      # 0.05 보다 작으면 등분산 가정 위배

예전 메모에는 stats.levene(residuals, y_pred) 로 적어 두었는데 이건 틀렸다. 레빈 검정은 여러 집단의 분산을 비교하는 것이지, 잔차와 적합값이라는 서로 다른 두 벡터를 넣는 검정이 아니다. 회귀의 등분산성은 Breusch-Pagan 이나 Goldfeld-Quandt 로 본다.

그림으로는 Scale-Location 플롯이다.

std_resid = model.get_influence().resid_studentized_internal
plt.scatter(model.fittedvalues, np.sqrt(np.abs(std_resid)))
plt.title('Scale-Location Plot')
plt.show()

점들이 깔때기 모양으로 벌어지면 이분산이다.

3. 독립성 — 더빈-왓슨

print(sm.stats.stattools.durbin_watson(residuals))

summary() 에도 이미 찍혀 나온다.

DW 값	뜻
2 에 가까움	자기상관 없음. 원하는 상태
0 ~ 2	양의 자기상관
2 ~ 4	음의 자기상관
0 또는 4 에 가까움	자기상관이 매우 강하다

시계열 데이터에 회귀를 적합하면 여기서 걸린다. DW 가 0 쪽으로 치우쳤다는 것은 회귀가 아니라 시계열 모형을 써야 한다는 신호다.

4. 정규성

sm.qqplot(residuals, stats.norm, fit=True, line='45')
print(stats.shapiro(residuals))     # p < 0.05 면 위배

정규성은 잔차에 대한 가정이지 원 데이터에 대한 가정이 아니다. $y$ 가 치우쳐 있어도 잔차가 정규면 문제없다.

5. 영향점 — Cook’s Distance

influence = model.get_influence()
(c, _) = influence.cooks_distance

outliers = np.where(c > 0.5)[0]      # 0.5 를 넘으면 영향점으로 본다
print("Outliers:", outliers)

fig, ax = plt.subplots(figsize=(12, 6))
ax.stem(c, linefmt="C0-", markerfmt="C0o", basefmt="k-")
ax.set_xlabel("Data points")
ax.set_ylabel("Cook's distance")
plt.show()

이상치와 영향점은 다르다. 혼자 멀리 떨어진 점(이상치)이라도 회귀선을 안 흔들면 놔둬도 되고, 값이 평범해 보여도 회귀선을 크게 돌리면(영향점) 문제가 된다. Cook’s D 는 뒤쪽을 잡는다.

변수 선택과 벌점화

21회 2번이 이 문제였다 — 선형회귀, 릿지, 라쏘를 각각 적합하고 최적 alpha 를 찾으라는.

from sklearn.linear_model import Ridge, Lasso
from sklearn.metrics import r2_score, mean_squared_error

best = (None, -np.inf)
for alpha in np.arange(0, 1.01, 0.1):
    m = Ridge(alpha=alpha).fit(X_train, y_train)
    score = m.score(X_test, y_test)
    if score > best[1]:
        best = (alpha, score)

print("최적 alpha:", best[0])
m = Ridge(alpha=best[0]).fit(X_train, y_train)
pred = m.predict(X_test)
print("R2:", r2_score(y_test, pred))
print("RMSE:", np.sqrt(mean_squared_error(y_test, pred)))

np.arange(0, 1.01, 0.1) 로 끝을 1.01 로 둔다. np.arange(0, 1, 0.1) 은 1.0 을 포함하지 않는다. “0부터 1까지 0.1 단위로 모두 탐색”이라는 문제 조건에서 마지막 값이 빠진다.

교차검증으로 한 번에 하려면 RidgeCV, LassoCV 다.

from sklearn.linear_model import LassoCV
from sklearn.preprocessing import StandardScaler
from sklearn.feature_selection import SelectFromModel

X_scaled = StandardScaler().fit_transform(X)

# 벌점화로 변수 선택 — 계수가 0 이 된 변수는 버린다
model = LassoCV(cv=5).fit(X_scaled, y)
coef = pd.Series(model.coef_, index=X.columns)
selected = coef[coef != 0].index

# 또는 중요도 기준으로 자른다
selection = SelectFromModel(LinearRegression().fit(X_scaled, y), threshold='median')
selection.fit(X_scaled, y)
selected = X.columns[selection.get_support()]

릿지와 라쏘 앞에는 반드시 표준화가 온다. 벌점이 계수의 크기에 걸리는데, 변수의 단위가 제각각이면 단위가 큰 변수만 벌을 받는다.

	벌점	성질
Ridge	$\lambda \sum \beta_j^2$	계수를 0 쪽으로 줄이지만 0 으로 만들진 않는다. 공선성 완화
Lasso	$\lambda \sum \lvert \beta_j \rvert$	계수를 정확히 0 으로 만든다. 변수 선택이 된다

Box-Cox 변환

종속변수가 치우쳐 있어 정규성이 깨질 때 쓴다. 17회의 log1p 도 같은 목적이다.

from scipy import stats

model = sm.OLS(y, sm.add_constant(x)).fit()
print(model.summary())

y_boxcox, lambda_ = stats.boxcox(y)          # y 는 모두 양수여야 한다
model_boxcox = sm.OLS(y_boxcox, sm.add_constant(x)).fit()
print(model_boxcox.summary())

변환 전후의 $R^2$ 와 잔차 QQ 플롯을 나란히 붙여 “좋아졌다”를 보인다. 변환 자체가 답이 아니라 변환의 효과를 보이는 것이 답이다. $\lambda$ 가 0 이면 로그 변환, 1 이면 변환 불필요와 같다는 것도 한 줄 써 준다.

AIC 와 BIC

$\text{AIC} = -2\ln L + 2k \qquad \text{BIC} = -2\ln L + k\ln n$

$-2\ln L$ 은 적합도, 뒤의 항은 파라미터 수 $k$ 에 대한 벌이다. 둘 다 작을수록 좋다.

변수를 늘리면 적합도는 항상 좋아진다. 그래서 적합도만으로 모형을 고르면 변수를 다 집어넣는 쪽이 이긴다. AIC·BIC 는 거기에 벌을 매겨서 적합도와 간결함의 균형점을 잡는다.

$n \ge 8$ 이면 $\ln n > 2$ 이므로 BIC 가 AIC 보다 변수 증가에 민감하다. 변수 개수를 줄이는 게 우선이면 BIC 를 본다.

OLS summary 읽는 법

print(model.summary()) 를 붙여 놓고 아래 항목을 문장으로 옮기면 그게 답안이다.

항목	뜻
R-squared / Adj. R-squared	설명력. 변수가 늘면 R²는 무조건 오르므로 Adj. 를 본다
F-statistic / Prob (F)	모형 전체의 유의성. $p < 0.05$ 면 “이 회귀식은 유의하다”
coef / P>\|t\|	각 계수와 그 유의성. $p \ge 0.05$ 인 변수는 기여가 없다
Omnibus / Prob(Omnibus)	잔차의 정규성. 왜도·첨도 기반. 값이 작을수록 좋다
Jarque-Bera (JB) / Prob(JB)	같은 목적의 다른 검정
Skew	왜도. 0 에 가까우면 대칭
Kurtosis	첨도. 3 에 가까우면 정규분포와 비슷. 크면 뾰족하다
Durbin-Watson	잔차의 자기상관. 2 근처면 독립
Cond. No.	다중공선성. 30 을 넘으면 의심한다

맨 아래 각주도 읽는다.

Standard Errors assume that the covariance matrix of the errors is correctly specified — 고전적 가정을 그대로 쓴 표준오차다
Standard Errors are robust to heteroscedasticity — 이분산이 있어도 견디는 방식으로 계산됐다
Standard Errors are clustered at the X level — X 기준으로 묶어서 계산했다

다항 회귀

21회 3번 — 3차까지 적합하고 차수별로 산점도·계수·회귀선을 그리라는 문제.

from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import make_pipeline

for degree in [1, 2, 3]:
    m = make_pipeline(PolynomialFeatures(degree), LinearRegression()).fit(X, y)
    xs = np.linspace(X.min(), X.max(), 200).reshape(-1, 1)
    plt.scatter(X, y)
    plt.plot(xs, m.predict(xs), color='red')
    plt.title(f'degree = {degree}')
    plt.show()
    print(degree, m.named_steps['linearregression'].coef_)

numpy 로 계수만 빨리 뽑을 수도 있다 — np.polyfit(x, y, 3).

베이지안 회귀

26회에 나왔다. 계수와 함께 예측의 불확실성을 준다는 게 보통 회귀와 다른 점이다.

from sklearn.linear_model import BayesianRidge

model = BayesianRidge().fit(X, Y)
print("alpha:", model.alpha_)          # 잡음의 정밀도
print("coef:",  model.coef_)
print("sigma:", np.sqrt(model.lambda_ / model.alpha_))

Y_pred, Y_std = model.predict(X_test, return_std=True)

plt.scatter(X, Y, label="Data")
plt.plot(X_test, Y_pred, color="red", label="Predicted Mean")
plt.fill_between(X_test.ravel(), Y_pred - 2*Y_std, Y_pred + 2*Y_std,
                 color="orange", alpha=0.3, label="Predicted Uncertainty")
plt.legend()
plt.show()

return_std=True 로 받은 표준편차로 ±2σ 띠를 칠하면 그림 하나로 설명이 끝난다.

분위 회귀

평균이 아니라 분위수를 적합한다. 이상치에 강하고, 조건부 분포의 모양을 볼 수 있다.

import statsmodels.api as sm

for q in [0.25, 0.5, 0.75]:
    model = sm.QuantReg(Y, sm.add_constant(X)).fit(q=q)
    plt.plot(X_pred, model.predict(sm.add_constant(X_pred)), label=f"Quantile {q}")

세 선이 나란하면 분산이 일정한 것이고, 벌어지면 이분산이다. 등분산성을 그림으로 보여주는 용도로도 쓸 만하다.