모형진단과 다항회귀
Contents
3.10. 모형진단과 다항회귀#
지금까지는 예측값이 직선형태로 나오는 선형(linear)회귀분석을 이용하여 종속변수를 예측하였다. 그러나 현실의 데이터는 이러한 선형회귀분석 모형을 따르지 않을 수도 있다. 이 절에서는 선형회귀분석에서 사용한 모형이나 가정이 적정한지 아닌지를 확인하고 그렇지 않은 이런 경우에 어떤 방식으로 모형을 수정해야하는지를 공부한다. 우리가 사용한 모형과 가정이 적정한지를 검토하는 작업을 모형진단(model diagnosis)이라고 한다.
3.10.1. 부분회귀#
독립변수가 하나뿐인 경우에는 x축을 독립변수로, y축을 종속변수로 하는 스캐터플롯과 예측값 직선을 이용하여 회귀분석의 결과를 시각적으로 표시하여 독립변수와 종속변수의 관계 및 회귀 성능을 유추할 수 있었다. 하지만 독립변수가 두개 이상인 경우에 이러한 방식을 사용하면 다른 독립변수의 영향으로 예측값이 직선으로 나오지 않기 때문에 이러한 시각화 방식을 사용하기 어렵다.
보스턴 집값 데이터의 일부를 lstat과 rm 두가지 독립변수로 선형회귀분석한 결과를 이용하여 예를 들어보자.
import statsmodels.api as sm
boston = sm.datasets.get_rdataset("Boston", "MASS").data.iloc[:100]
model = sm.OLS.from_formula("medv ~ scale(lstat) + scale(rm)", boston)
result = model.fit()
print(result.summary())
OLS Regression Results
==============================================================================
Dep. Variable: medv R-squared: 0.799
Model: OLS Adj. R-squared: 0.795
Method: Least Squares F-statistic: 193.1
Date: Sat, 06 Aug 2022 Prob (F-statistic): 1.52e-34
Time: 16:05:07 Log-Likelihood: -239.08
No. Observations: 100 AIC: 484.2
Df Residuals: 97 BIC: 492.0
Df Model: 2
Covariance Type: nonrobust
================================================================================
coef std err t P>|t| [0.025 0.975]
--------------------------------------------------------------------------------
Intercept 22.3090 0.268 83.135 0.000 21.776 22.842
scale(lstat) -2.1390 0.339 -6.319 0.000 -2.811 -1.467
scale(rm) 3.6891 0.339 10.898 0.000 3.017 4.361
==============================================================================
Omnibus: 5.327 Durbin-Watson: 1.248
Prob(Omnibus): 0.070 Jarque-Bera (JB): 4.802
Skew: 0.419 Prob(JB): 0.0906
Kurtosis: 3.670 Cond. No. 2.03
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
예를 들어 x축을 첫번째 독립변수인 lstat으로, y축을 종속변수인 medv로 놓고 2차원 평면에 데이터(o)와 예측값(x)을 그리면 다음과 같다.
x = boston.lstat
y = boston.medv
y_hat = result.fittedvalues
plt.figure(figsize=(10, 8))
plt.scatter(x, y, marker="o", label="원 데이터")
plt.scatter(x, y_hat, marker="x", label="예측값")
plt.legend()
plt.title("다중회귀분석 결과를 2차원에 표시한 경우")
plt.show()
statsmodels 패키지가 제공하는 plot_fit 함수를 사용하면 좀 더 쉽게 그릴 수 있다. 이 때는 각 데이터에 대한 개별 잔차의 표준편차를 바(bar)로 표시해 준다. 개별 잔차의 표준편차보다도 바깥에 있는 데이터는 아웃라이어일 가능성이 높다.
f, ax = plt.subplots(figsize=(10, 8))
sm.graphics.plot_fit(result, "scale(lstat)", ax=ax)
plt.show()
이러한 경우에 사용할 수 있는 방법이 CCPR(component-plus-residual) 및 부분회귀(partial-regression) 시각화 방법이다.
CCPR(component-plus-residual) 시각화는 2차원 평면에서 x축을 특정한 독립변수 하나, 예를 들어 \(k\)번째 독립변수인 \(x_k\) 값으로 놓고 y축은 다음 데이터를 이용한다.
원 데이터의 y축 값으로는
를 사용한다. 여기에서 \(e\)는 전체 독립변수를 사용한 다중선형회귀분석 결과에서 나온 잔차다.
예측 데이터의 y축 값으로는 선형회귀 예측값에서 \(x_k\)과 관련된 부분인
만 사용한다. 이렇게 하면 다중선형회귀분석의 결과로 나온 가중치의 값(기울기)도 표시하고 그에 따른 오차도 표시된다.
statsmodels 패키지는 CCPR 시각화를 위한 plot_ccpr 함수를 제공한다. 사용법은 다음과 같다.
f, ax = plt.subplots(figsize=(10, 8))
sm.graphics.plot_ccpr(result, "scale(lstat)", ax=ax)
plt.show()
또다른 방법은 부분회귀 방법을 사용하는 것이다. 부분회귀(partial regression)는 특정한 독립변수 \(x_k\)과 종속변수 \(y\)만의 관계만 분석하고 \(x_k\)를 제외한 다른 독립변수의 영향은 제거하는 방법이다. 부분회귀를 하려면 원리적으로 다음과 같이 3번의 회귀분석을 해야 한다.
특정한 독립변수 \(x_k\)를 제외한 나머지 독립변수들로 종속변수 \(y\)를 선형회귀분석하여 잔차 \(y_k^{\ast}\)를 구한다.
특정한 독립변수 \(x_k\)를 제외한 나머지 독립변수들로 특정한 독립변수 \(x_k\)를 선형회귀분석하여 잔차 \(x_k^{\ast}\)를 구한다.
잔차 \(x_k^{\ast}\)를 독립변수로, 잔차 \(y_k^{\ast}\)를 종속변수로 하여 선형회귀분석한다.
예를 들어 첫번째 독립변수인 \(x_1\)과 종속변수 \(y\)만의 관계만 구하고 싶은 경우, 1번 작업을 통해 \(x_1\)을 제외한 \(x_2, x_3, \ldots, x_k\)만 사용하여 선형회귀분석을 한다. 그러면 이 때 나오는 잔차 \(y_k^{\ast}\)는 독립변수 \(x_2, x_3, \ldots, x_k\)로는 도저히 예측할 수 없는 부분만 남는다. 그리고 2번 작업을 통해 \(x_1\) 정보 중 \(x_2, x_3, \ldots, x_k\)와 상관관계가 있는 값을 제거하여 순수한 \(x_1\) 정보만을 남긴다. 마지막으로 3번 작업을 통해 순수한 \(x_1\) 정보와 \(x_2, x_3, \ldots, x_k\)로 예측 불가능한 정보 사이의 관계를 구하는 것이다. 이러한 방식으로 부분회귀를 하면 3번 작업에서 나온 가중치 \(\hat{w}_k\)가 모든 독립변수를 사용하였을 때 나오는 가중치 \(\hat{w}_k\)와 일치함을 수학적으로 증명할 수 있다.
statsmodels 패키지는 부분회귀 시각화를 위한 plot_partregress 함수를 제공한다. 이 함수를 사용할 때는 다음과 같은 인수를 넣어야 한다.
종속변수의 이름 문자열
분석하고자 하는 특정 독립변수의 이름 문자열
분석하고자 하는 특정 독립변수를 제외한 나머지 독립변수들의 이름 문자열 리스트
데이터프레임
obs_labels: 각 데이터에 라벨을 그리는 경우True, 아니면False
함수를 실행하면 x축에는 잔차 \(x_k^{\ast}\)를, y축에는 잔차 \(y_k^{\ast}\)를 스케터플롯으로 그리고 가중치 \(\hat{w}_k\)를 사용한 선형회귀 결과를 직선으로 표시한다.
f, ax = plt.subplots(figsize=(10, 8))
sm.graphics.plot_partregress("medv", "scale(lstat)", ["scale(rm)"], boston,
obs_labels=False, ax=ax)
plt.show()
eval_env: 1
plot_regress_exog 함수를 이용하면 CCPR과 부분회귀 결과를 동시에 그릴 수 있다.
f, ax = plt.subplots(figsize=(10, 8))
sm.graphics.plot_regress_exog(result, "scale(lstat)", fig=f)
plt.tight_layout()
plt.show()
eval_env: 1
상단 왼쪽의 그림은 우리가 처음에 그렸던 plot_fit 함수의 결과다. 하단은 부분회귀와 CCPR 결과다.
plot_regress_exog 함수가 출력한 그림 중 상단 오른쪽의 그림은 잔차플롯(residual plot)이라고 하여 독립변수와 그에 대응하는 잔차의 값을 그린 것이다.
3.10.2. 잔차분석#
잔차를 분석하는 이유는 우리가 만든 선형회귀분석 모형이 정상적인 모형인지 아닌지를 검증하기 위한 것이다. 만약 우리가 사용한 모형과 데이터가 선형회귀분석에 사용되는 기본적인 확률론적 가정을 만족한다면 잔차는 평균이 0이고 분산이 고정된 정규분포를 따라야 한다.
회귀분석 결과 요약에서 하단에 있는 부분이 잔차의 정규분포에 대한 통계치 및 정규성 검정 결과다.
Omnibus: 5.327
Prob(Omnibus): 0.070 Jarque-Bera (JB): 4.802
Skew: 0.419 Prob(JB): 0.0906
Kurtosis: 3.670
각 항목의 의미는 다음과 같다.
Skew: 잔차분포의 왜도(skewness)Kurtosis: 잔차분포의 첨도(kurtosis)Omnibus: 옴니버스 정규성 검정의 통계치Prob(Omnibus): 옴니버스 정규성 검정의 유의확률Jarque-Bera (JB): 자끄베라 정규성 검정의 통계치Prob(JB): 자끄베라 정규성 검정의 유의확률
우리가 사용한 예제에서는 정규성 검정결과가 다음과 같다.
옴니버스 정규성 검정의 유의확률 = 7%
자끄베라 정규성 검정의 유의확률 = 9%
만약 유의수준이 5%라면 유의확률이 이보다 크므로 정규성 검정을 통과한다. 즉 모형의 기본가정을 만족한다.
만약 잔차가 정규성을 만족하지 못하면 우리가 지금까지 사용한 가중치의 유의확률이나 신뢰구간 등의 분석결과에 오차가 있을 수 있다. 따라서 원칙적으로는 정규성 분포를 만족할 때까지 아웃라이어를 제거하거나 모형을 수정할 필요가 있다.
그러나 현실적으로 데이터의 개수가 충분한 경우(보통 15개 이상), 잔차의 비정규성이 분석결과에 큰 영향을 미치는 경우는 드물고 대부분의 문제에서는 비실용적일 정도로 복잡한 모형이 아닌 한 정규성 검정읕 통과하지 못하는 경우가 많다. 정규성 검정의 중요성에 대해서는 아직도 통계학자들 사이에서 논란이 되고 있는 부분이지만 이 책에서는 일단 이 부분을 크게 신경쓰지 않고 넘어가기로 한다.
선형회귀에서 실질적으로 정규성보다 더 중요한 잔차의 성질은 다음 성질이다.
잔차의 조건부 평균
선형회귀분석으로 나온 \(N\)개의 잔차 \(e\)의 기댓값은 독립변수와 상관없이 항상 0이다.
잔차의 조건부 분산
선형회귀분석으로 나온 \(N\)개의 잔차 \(e\)의 분산은 독립변수와 상관없이 항상 일정하다. 이를 등분산성(homoscedasticity)이라고 한다. 만약 이를 만족하지 않으면 이분산성(heteroscedasticity)이라고 한다.
statsmodels 패키지에는 잔차의 조건부 평균이 0인지 검증하기 위한 통계정 검정 기능이 없다 따라서 독립변수와 잔차의 관계를 나타내는 잔차플롯을 그려 시각적으로 살펴볼 수 밖에 없다. 안타깝게도 statsmodels 패키지는 잔차플롯만 그리는 함수도 없다. 따라서 잔차플롯만 그리고 싶을 때는 저자가 만든 다음 코드를 사용한다. 이 잔차플롯에서 파란색 곡선은 독립변수의 변화에 따른 잔차의 평균값을 표시한 것이다.
def resid_plot(result, exo_name):
id_x = result.model.exog_names.index(exo_name)
x = result.model.exog[:, id_x]
r = result.resid
l = sm.nonparametric.lowess(r, x, frac=0.3)
f, ax = plt.subplots(figsize=(10, 8))
plt.scatter(x, r, c="r", alpha=0.5)
plt.plot(l[:, 0], l[:, 1], lw=5, alpha=0.8)
plt.axhline(y=0, color='black')
plt.ylim(-10, 10)
plt.title(f"잔차플롯: 독립변수={exo_name}")
resid_plot(result, "scale(lstat)")
plt.show()
잔차 평균이 독립변수와 관계없이 0이 되려면 잔차플롯에서 푸른 곡선이 y=0 에 해당하는 평평한 직선이어야 한다. 그러나 독립변수와 잔차의 관계를 보면 독립변수 값이 너무 작거나 커지면 잔차의 값이 양의 방향으로 증가하는 것을 볼 수 있다. rm 독립변수에 대해 잔차플롯을 그려도 비슷한 현상이 발생한다.
resid_plot(result, "scale(rm)")
잔차의 조건부 분산 성질의 경우에는 statmodels 패키지에 이를 검증하기 위한 다음과 같은 함수가 있다.
het_breuschpagan: 이분산성에 대한 브레슈-페이건(Breusch-Pagan) 검정het_white: 이분산성에 대한 화이트(white) 검정
이 검정들의 귀무가설은 등분산성이므로 유의확률이 큰 경우에는 등분산성이고 유의확률이 작으면 이분산성이다.
from statsmodels.stats.diagnostic import het_breuschpagan, het_white
이 함수들의 반환값은 다음과 같은 순서로 나오는 값의 튜플이다.
라그랑주 승수검정 통계치
라그랑주 승수검정 유의확률
F검정 통계치
F검정 유의확률
두번째와 4번째의 유의확률의 값이 크면 등분산성이다. 우리가 사용한 예제에서는 이분산성이 존재하는 것을 확인할 수 있다.
het_breuschpagan(result.resid, model.exog)
(22.398766185277374,
1.3682634380455558e-05,
13.99900628615844,
4.55646855722389e-06)
het_white(result.resid, model.exog)
(23.458466930318266,
0.00027580773458726204,
5.7618283904568415,
0.00011057945959355164)
이분산성이 존재하는 경우 가중치의 유의확률이나 신뢰구간에 오차가 있을 수 있으므로 앞에서 말한바와 마찬가지로 아웃라이어를 제거하거나 모형을 수정하여 다시 검정을 실시할 필요가 있다.
3.10.3. 다항회귀#
잔차분석을 했을 때 잔차가 원래 가져야할 성질을 만족하지 않는 이유는 우리가 사용한 회귀모형이 데이터에 맞지 않기 때문이다. 즉 데이터가 다음과 같은 선형관계를 따르지 않는다
따라서 선형모형을 확장한 비선형 모형을 사용해야한다. 그런데 일반적인 비선형모형을 사용하면 지금까지 우리가 공부한 선형회귀분석을 활용할 수 없기 때문에 다음과 같은 비선형모형을 사용하는 것이 좋다.
이 식에서 \(f_l()\)은 독립변수 \(x_1, x_2, \ldots, x_K\)를 입력으로 가지는 비선형함수들이다.
이 모형은 다음과 같이 표시할 수도 있다.
즉 기존에 우리가 사용한 독립변수 \(x_1, x_2, \ldots, x_K\)를 비선형 변환한 값을 새로운 독립변수 \(x'_1, x'_2, \ldots, x'_L\)로 사용하는 것이다. 이 때 비선형함수는 자유롭게 지정할 수 있기 때문에 새로운 독립변수의 개수 \(L\)는 기존의 독립변수의 개수 \(K\)보다 훨씬 많아질 수 있다.
이 방식을 사용하면 실질적으로는 비선형 함수를 사용하면서 우리가 지금까지 공부한 선형회귀분석을 활용할 수 있다.
이 때 비선형함수 \(f_l()\)를 만드는 방법은 데이터 분석가의 마음대로지만 가장 간단한 방법은 개별독립변수의 고차항 \(x^2, x^3, x^4, \ldots\)를 사용하는 방법이다. 예를 들어 독립변수가 하나뿐인 경우 회귀모형은 다음과 같아진다.
이 방법을 다항회귀(polynomial regression)이라고 한다.
statsmodels 패키지에서 다항회귀를 하려면 모형 문자열에 거듭제곱을 나타내는 ** 연산과 I() 연산을 사용하여 독립변수를 추가한다.
formula = "y ~ x + I(x**2) + I(x**3) + I(x**4) + I(x**5)"
만약 I() 연산을 하지 않으면 나중에 설명할 상호작용(interaction) 항이 추가되어 다른 모형이 되므로 주의해야 한다.
예제로 사용한 보스턴 집값 분석에서 rm과 lstat의 2차항을 추가하여 분석하면 다음과 같다.
formula2 = "medv ~ scale(lstat) + scale(rm) + scale(I(lstat**2)) + scale(I(rm**2))"
model2 = sm.OLS.from_formula(formula2, boston)
result2 = model2.fit()
print(result2.summary())
OLS Regression Results
==============================================================================
Dep. Variable: medv R-squared: 0.853
Model: OLS Adj. R-squared: 0.847
Method: Least Squares F-statistic: 137.8
Date: Sat, 06 Aug 2022 Prob (F-statistic): 1.15e-38
Time: 16:05:09 Log-Likelihood: -223.48
No. Observations: 100 AIC: 457.0
Df Residuals: 95 BIC: 470.0
Df Model: 4
Covariance Type: nonrobust
========================================================================================
coef std err t P>|t| [0.025 0.975]
----------------------------------------------------------------------------------------
Intercept 22.3090 0.232 96.161 0.000 21.848 22.770
scale(lstat) -6.1843 0.998 -6.195 0.000 -8.166 -4.203
scale(rm) -14.1352 4.350 -3.250 0.002 -22.770 -5.500
scale(I(lstat ** 2)) 3.4907 0.915 3.813 0.000 1.673 5.308
scale(I(rm ** 2)) 17.0996 4.304 3.973 0.000 8.555 25.645
==============================================================================
Omnibus: 25.903 Durbin-Watson: 1.424
Prob(Omnibus): 0.000 Jarque-Bera (JB): 65.023
Skew: 0.899 Prob(JB): 7.59e-15
Kurtosis: 6.517 Cond. No. 46.2
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
일단 가중치의 유의확률로 볼 때 모든 독립변수가 유의하고 결정계수 성능은 증가한 것을 알 수 있다.
잔차플롯을 그려보면 잔차의 평균이 앞의 모형보다 0에 더 가까운 형태로 개선되었다.
resid_plot(result2, "scale(lstat)")
resid_plot(result2, "scale(rm)")
resid_plot(result2, "scale(I(lstat ** 2))")
resid_plot(result2, "scale(I(rm ** 2))")
이분산성 검정을 한 결과도 유의확률이 1% 이상으로 증가하여 유의수준이 1%이면 검정을 통과하여 등분산성을 보이게 된다.
het_breuschpagan(result2.resid, model2.exog)
(10.32306667382783,
0.0353234696975647,
2.7339564859078136,
0.03337658699888184)
het_white(result2.resid, model2.exog)
(13.899283201672985,
0.30718609379310136,
1.1703712461321483,
0.3170099161339497)