다중 회귀모형의 진단
내용
단순회귀모형과 같이 회귀분석은 몇 가지 가정을 전제한 상태에서 모형을 구축합니다. 그러므로 그 가정의 충족되지 않은 자료들에서 생성된 모형의 경우 적용과 적합성에 문제가 발생됩니다. 예를 들어 데이터가 정규분포를 따르지 않는 상태에서 구축된 회귀모형의 예측 범위는 매우 넓어질 수 있으므로 그 자체의 의미가 감소됩니다. 또한 변수들의 독립성에 문제가 있는 경우 반응변수에 대한 설명변수들의 선택에서 발생하는 오류는 실제와 다른 결과를 가져올 수 있습니다. 이 원인들은 모두 회귀모형의 신뢰도를 악화시키며 예측정도를 빈약하게 만들 수 있습니다.
정규성(normality) 과 이상치(outlier) 평가
예 1)
다음 코드에 의해 생성된 자료 코스피지수(kos), 코스탁지수(kq), kodex 레버리지(kl), kodex 인버스(ki), 그리고 원달러환율(WonDol)의 일일 종가들을 설명변수로 사용하여 삼성전자(sam)의 일일 종가를 추정하는 회귀모델을 구축해 봅니다. (이 자료에서 설명변수는 반응변수보다 1일 앞선 데이터 입니다.)
import numpy as np
import pandas as pd
import yfinance as yf
st=pd.Timestamp(2023,1, 10)
et=pd.Timestamp(2024, 5, 30)
code=["^KS11", "^KQ11", "122630.KS", "114800.KS","KRW=X","005930.KS"]
nme=["kos","kq","kl", "ki", "WonDol","sam" ]
da={}
for i, j in zip(nme,code):
da[i]=yf.download(j,st, et)["Close"]
da2=pd.DataFrame(da.values(), index=da.keys()).T
da2=da2.ffill()
da2.head(3).round(0)
자료를 표준화하여 ols() 클래스를 사용 중회귀모형을 생성합니다(이 모형의 생성과정은 다중회귀모델의 생성 참조).
from sklearn import preprocessing import matplotlib.pyplot as plt from statsmodels.formula.api import ols from scipy import stats
ind=da2.values[:-1,:-1] de=da2.values[1:,-1].reshape(-1,1) final=da2.values[-1, :-1].reshape(1,-1) indScaler=preprocessing.StandardScaler().fit(ind) deScaler=preprocessing.StandardScaler().fit(de) indNor=indScaler.transform(ind) finalNor=indScaler.transform(final) deNor=deScaler.transform(de) da3=pd.DataFrame(np.c_[indNor, deNor.reshape(-1,1)]) da3.columns=da2.columns form='sam~kos+kq+kl+ki+WonDol' reg=ols(form, data=da3).fit()
반응 변수가 고정된 설명 변수 집합에 대해 정규 분포를 따르는 경우 잔차 값은 평균이 0인 정규 분포를 따라야 합니다.
자료의 정규성은 모델에 의해 생성되는 잔차의 분석으로 검정할 수 있습니다. 잔차와 스튜던트 잔차는 OLS 모델의 .resid, .residpearson 속성으로 확인할 수 있습니다. 잔차의 정규성은 shapiro()를 사용합니다.
err=reg.resid
normT=stats.shapiro(err)
print("f-static: %.3f, p-value: %.3f" %(normT[0], normT[1]))
f-static: 0.943, p-value: 0.000
위 결과에 의하면 잔차는 정규분포에 부합한다는 귀무가설을 기각할 수 있습니다. 그림 1은 오차에 대한 Q-Q plot로서 양 끝단 값들에 의한 기울기의 현저한 이탈정도는 위 결과를 부연합니다.
plt.figure(figsize=(4,3)) stats.probplot(err, plot=plt, rvalue=True) plt.show()
모델의 잔차는 위 결과와 같이 정규성에 부합합니다. 이것은 그 데이터에 이상치에 의한 영향이 매우 작다는 것을 의미합니다. 이상치를 결정하기 위해 스튜던트 잔차(rstudent)와 Cook's Distance(D)를 확인합니다. 그림 2는 이상치를 판단하는 기준인 스튜던트 표준잔차(rstudent)를 나타낸 것입니다.
fig=plt.figure(figsize=(4, 3))
plt.stem(reg.resid_pearson)
plt.axhline(3, c="g", ls="--")
plt.axhline(-3, c="g", ls="--")
plt.title("studented Residuals")
plt.show()
스튜던트 잔차의 경우 [-3, 3] 구간외에 존재하는 데이터의 경우 이상치로 고려합니다. 위 경우는 이상치로 의심할 데이터 존재를 나타냅니다.
hat 행렬을 기반으로 하는 레버리지, cook's distance(CD)를 사용하여 이상치를 결정하기 위한 계산은 OLS 모델의 get_influence() 메소드로 계산됩니다.
inf=reg.get_influence().summary_frame() inf["cooks_d"].head(3)
0 0.002426 1 0.004073 2 0.003554 Name: cooks_d, dtype: float64
cd_idx=np.where(inf["cooks_d"]>0.5)[0] print(cd_idx)
[]
일반적으로 CD값이 0.5 이상이면 이상치일 가능성을 보입니다. 이 자료의 경우 이상치는 없는 것으로 판단할 수 있습니다. 그러나 스튜던트 잔차에 의한 잠재적 이상치를 제거한 후 회귀모델을 구축하고 오차의 정규성을 확인해봅니다.
resP_id=np.where(reg.resid_pearson>3)[0] print(resP_id)
[130 131 132]
ind_resP=np.delete(indNor, resP_id, axis=0)
de_resP=np.delete(deNor, resP_id, axis=0)
reg_resP=sm.OLS(de_resP, ind_resP).fit()
print(f"회귀계수: {np.around(reg_resP.params, 3)}\nR2: {reg_resP.rsquared.round(3)}")
회귀계수: [-2.51 0.449 1.142 -1.873 0.076] R2: 0.893
stat, pVal=stats.shapiro(reg_resP.resid) np.around(pd.DataFrame([stat,pVal], index=["통계량", "p-value"]), 3)
| 0 | |
|---|---|
| 통계량 | 0.988 |
| p-value | 0.004 |
plt.figure(figsize=(3,2)) fig=stats.probplot(reg_resP.resid, plot=plt) plt.show()
위 결과에 의하면 이상치제거에 의해 약간의 개선 효과가 일어남을 확인할 수 있습니다. 그러나 정량적으로 정규성에 부합한다는 귀무가설을 기각할 수 있습니다.
독립성
설명변수들은 서로 독립이어야 합니다. 그러나 이것에 대한 정량적 평가는 모호하며 데이터 자체와 수집 등에 대한 이해에 기반한 정성평가가 우선됩니다. 그러나 시간에 따른 데이트들의 상관정도는 정량화할 수 있습니다.
선형성과 자기상관
오차는모델에 의한 예측지와 실측치의 차이이므로 회귀계수와 확률변수들에 의해 결정됩니다. 그러므로 오차 역시 확률분포로서 어떤 체계적인 경향을 발견될 수 없습니다. 그러나 모형에 참여하는 각 변수에서 자기상관을 보이는 경우 오차 역시 특정한 경향을 보일 수 있습니다. 특히 금융자료와 같이 시계열 자료(시간에 따른 자료)의 경우 자기상관 문제가 두드러집니다. 즉, 시간적으로 더 가깝게 수집된 관찰은 시간적으로 먼 관찰보다 서로 더 상관관계가 있습니다. 자기상관 검정은 Durbin-Watson 테스트의 결과로 판단할 수 있습니다. 이 결과는 [0, 4] 범위의 값을 반환하며 2를 기준으로 낮으면 양의 자기상관 반대는 음의 자기 상관을 나타냅니다.
sm.stats.stattools.durbin_watson(err).round(3)
0.625
위 결과는 (1.5, 2.5) 구간을 벗어난 값으로 양의 자기상관을 의심할 수 있습니다. 자기상관성을 가진 데이터의 경우 설명변수에 일정한 가중치를 적용하는 lasso, ridge 등을 적용합니다.
등분산성
오차의 등분산성을 검정하기 위한 방법인 Breusch-Pegan 검정을 소개합니다. 이 검정의 귀무가설(H0)는 "등분산성의 성립"입니다.
Breusch-Pegan 검정의 기본 아이디어는 설명변수를 기준으로 두개 이상의 그룹으로 구분하고 각 그룹의 잔차 제곱의 합 즉 SSE에 대해 모든 그룹의 합은 χ2분포를 따르므로 이 분포를 근거로 검정합니다. χ2분포에 부합하면 각 그룹은 정규분포를 따르고 독립적이므로 각 분포의 차이가 없다는 것을 의미합니다. 즉, 등분산성이 성립함을 나타냅니다. 그러므로 Breusch-Pegan 검정의 귀무가설(H0)은 '등분산성의 성립'이 됩니다. 이 검정은 다음과정으로 실행합니다.
- 설명변수와 반응변수 사이의 회귀모델 생성: regModel = 반응변수 ~ 설명변수
- regModel로 부터의 발생되는 잔차의 제곱에 대해 회귀모델 생성 : regModel2 = regModel.resid2 ~ 설명변수
- regModel2.resid에 대해 χ2 검정을 실시
위 과정의 3번째 단계는 statsmodels.stats.diagnostic.het_breuschpagan(잔차, 설명변수) 함수를 적용하여 검정 결과를 확인할 수 있습니다.
import statsmodels.api as sm
indNor0=sm.add_constant(indNor)
reg_res=sm.OLS(err**2, indNor0).fit()
lmv, lmp, fv, fp=sm.stats.diagnostic.het_breuschpagan(reg_res.resid,reg_res.model.exog)
print('라그랑주승수통계량: %.3f, p-value: %.3f, f-value: %.3f, fp-value: %.3f' %(lmv, lmp, fv, fp))
라그랑주승수통계량: 26.255, p-value: 0.000, f-value: 5.569, fp-value: 0.000
위 결과는 모델의 잔차는 등분산성을 따른다는 귀무가설을 기각할 수 있습니다. 그러나 스튜던트 잔차에 의한 잠재 이상치를 제거한 후의 모델에 대한 잔차의 등분산성은 개선 효과가 나타나며 유의 수준 0.05 기준으로 등분산성에 대한 귀무가설을 채택할 수 있습니다.
indNor1=sm.add_constant(ind_resP)
reg_res2=sm.OLS(reg_resP.resid**2, indNor1).fit()
lmv, lmp, fv, fp=sm.stats.diagnostic.het_breuschpagan(reg_res2.resid,reg_res2.model.exog)
print('라그랑주승수통계량: %.3f, p-value: %.3f, f-value: %.3f, fp-value: %.3f' %(lmv, lmp, fv, fp))
라그랑주승수통계량: 10.877, p-value: 0.054, f-value: 2.206, fp-value: 0.053