다중 공선성(Multicolinearity)
최소자승법은 모델에 의해 생성되는 오차를 최소화하도록 설계된 방법입니다. 식 1은 기사 '회귀계수의 추정: 최소제곱법(Least Square method)'의 식 8을 자세히 나타낸 것입니다.
\begin{align} \text{MSE}&=(y-X\beta)^T(y-X\beta)\\ \frac{\partial \text{MSE}}{\partial \beta}&=\frac{\partial }{\partial \beta}(y-X\beta)^T(y-X\beta)=0\\ & \Leftrightarrow \frac{\partial }{\partial \beta}(y^T-X^T\beta^T)(y-X\beta)\\ \tag{식 1}& \Leftrightarrow \frac{\partial }{\partial \beta}\left(y^Ty-y^TX\beta -\beta^TX^Ty+\beta^TX^TX\beta \right) \\ & \Leftrightarrow \frac{\partial }{\partial \beta}\left(y^Ty-y^TX\beta -(y^TX\beta)^T+\beta^TX^TX\beta \right) \\ & \Leftrightarrow -y^TX - X^Ty + 2X^TX\beta =0\\ & \Leftrightarrow X^TX\beta = X^Ty \\ & \Leftrightarrow \beta=(X^TX)^{-1}X^Ty\\ \because&\; X^Ty \leftrightarrow y^TX,\quad b^TX^TXb \leftrightarrow X^2b^2\end{align}
식 1에서 나타낸 것과 같이 회귀계수는 XTX의 역행렬에 의존합니다. 또한 행렬 XTX에 의해 공분산 행렬를 계산할 수 있습니다(공분산과 상관계수의 식 6 참조). 그러므로 이 행렬의 대각요소들을 각 변수의 분산, 대각외요소들을 공분산의 시작으로 고려할 수 있습니다. 행렬 XTX가 나타내는 분산과 공분산 지표와 역행렬 생성과정을 연결하여 봅시다.
식 2는 간단한 2×2 행렬의 역행렬을 구하는 방식을 나타내며 이 식으로부터 역행렬의 크기는 우항의 첫번째 항인 행렬식에 의존함을 알 수 있습니다. 즉, 행렬의 대각요소들의 곱과 대각외요소들의 곱의 차이가 클수록 역행렬의 요소들의 크기가 작아집니다. 대각외요소들은 공분산을 의미하는 지표가 되므로 공분산이 작아지면 역행렬의 크기가 작아집니다. 공분산이 작아진다는 것은 상관계수가 작다는 것을 의미합니다.
$$\tag{식 2} \begin{bmatrix}a& b \\ c& d \end{bmatrix}^{-1} = \frac{1}{ad-bc}\begin{bmatrix}d& -b \\ -c & a \end{bmatrix}$$
행렬 (XTX)-1와 오차의 분산을 사용하여 회귀계수의 분산을 계산할 수 있습니다. 위에서 언급한 것과 같이 (XTX)-1의 규모는 변수들간의 상관계수에 의존합니다. 큰 상관계수일 경우 역행렬의 크기가 증가하므로 회귀계수의 분산 역시 증가합니다. 이것은 추정의 신뢰구간을 증가시킵니다. 이와 같이 각 변수들 사이의 상관성이 증가함으로서 모델의 정확도를 감소시키는 현상을 다중공선성이라 합니다.
예 1)
다음은 설명변수들의 상관계수가 낮은 경우, 1인 경우 그리고 높은 경우에 따라 회귀계수 분산을 계산한 것입니다.
다음 작성한 함수 cofVar()은 설명변수(X)와 반응변수(y)에 대한 회귀모델을 생성하고 그 모델의 회귀계수의 분산을 반환합니다. 이 예제의 3가지 경우에 대한 회귀계수의 분산을 확인하기 위해 사용합니다.
def cofVar(X, y):
n, p=X.shape
re={}
m=LinearRegression().fit(X, y)
res=y-m.predict(X)
mse=np.sum(res**2)/(n-p-1)
re["cor"]=np.corrcoef(X, rowvar=False)
re["coef"]=m.coef_
re["mse"]=mse
try:
re["coef_var"]=la.inv(X.T@X)*mse
except:
re["coef_var"]=la.pinv(X.T@X)*mse
return(re)
X=stats.norm.rvs(size=[10, 2], random_state=1) X1=np.c_[X[:,0], 2*X[:,0]] X2=np.c_[X[:,0], X[:,0]+np.random.rand(10)] y=1.5*np.mean(X, axis=1)+np.random.rand(10) m1=cofVar(X,y) m2=cofVar(X1,y) m3=cofVar(X2, y)
cor=np.array([m1["cor"][0,1],m2["cor"][0,1],m3["cor"][0,1]]) mse=np.array([m1['mse'],m2['mse'],m3['mse']]) coef_var=np.array([m1['coef_var'][1,1],m2['coef_var'][1,1],m3['coef_var'][1,1]]) np.around(pd.DataFrame(np.c_[cor,mse,coef_var], columns=["cor","mse","var_b"]),3)
| cor | mse | var_b | |
|---|---|---|---|
| 0 | -0.361 | 0.102 | 0.012 |
| 1 | 1.000 | 0.571 | 0.009 |
| 2 | 0.950 | 0.549 | 0.181 |
위 결과와 같이 변수들 사이의 상관계수가 증가할수록 잔차의 분산(MSE)이 증가합니다. 위 결과의 두 번째 상관계수가 1인 경우는 설명변수와 반응변수가 정확히 선형관계 있는 것으로 계수의 분산은 낮을 수 밖에 없습니다. 이러한 특수한 경우를 제외하고 일반적으로는 위의 첫번째와 세번째 행의 결과와 같이 상관성의 증가는 MSE와 계수 분산의 증가를 유발하므로 각 계수의 신뢰구간이 넓어집니다. 결과적으로 모델에 의한 예측 구간이 넓어지므로 정확도 저하의 원인이 됩니다. 특히 이러한 낮은 정확도는 모델 구축에 사용되지 않은 새로운 데이터(검정데이터)에 대해서 두드러지게 나타날 수 있습니다. 즉, 과적합(overfitting)의 주요 원인이 됩니다.
모델을 생성하는 경우 자료를 모델 생성에 참여하는 세트와 참여하지 않고 생성된 모델을 검정하는 과정에서만 사용되는 세트로 구분합니다. 전자를 훈련데이터(training dataset), 후자를 검정데이터(test dataset)라고 합니다. 사용할 데이터가 풍부할 경우에 사용하며 확률·통계의 전제조건에서 비교적 자유로운 기계학습(machine leaning)등에서는 분석의 필수적인 단계입니다. 훈련데이터에서 보이는 정확도에 비해 검정데이터에 낮은 정확도를 과적합(overfitting)이라 합니다.
높은 상관성을 가지는 변수들 중의 하나만 선택하는 경우 다중공선성으로 인한 문제는 감소됩니다. 예를 들면 높은 상관성을 가지는 두 변수 중 하나의 변수만을 적용할 경우 모델의 적합도(R2)는 변화는 작을 것입니다. 이와 같이 특정한 변수를 제외할 경우 R2의 변화를 관찰하는 것으로 변수를 선택할 수 있습니다. 이러한 아이디어로 사용되는 결정기준을 분사팽창인수(varinacne inflation factor, VIF)라 하며 식 3과 같이 정의됩니다.
$$\tag{식 3}\text{VIF}_i=\frac{1}{1-R_i^2}$$
위 식에서 R2i는 i번째 변수를 제외한 나머지를 설명변수, i 번째 변수를 반응변수로 회귀분석시 결정계수입니다. 예로 VIF = 10인 경우 R2 = 0.9가 됩니다. 즉, i 번째 변수를 제외해도 모델은 반응변수를 90% 설명할 수 있음을 의미합니다. 다시 말하자면 i 번째 변수는 모델의 설명력에 미치는 영향이 작기 때문에 제외할 수 있음을 의미합니다.
일반적으로 다중 공선성에 대한 VIF의 기준을 표 1과 같습니다.
| 범위 | 상관성 |
|---|---|
| 1 = vif | 없음 |
| 1 < vif <5 | 낮음 |
| 5 ≤ vif | 높음 |
VIF는 statsmodel.stats.outliers_influence.variance_inflation_factor(변수, 제외할 변수의 인덱스)함수를 적용하여 계산할 수 있습니다. 또한 sklearn.LinearRegression() 클래스, 모델.score() 함수들을 적용하여 계산 할 수 있습니다.
예 2)
6개 주식에 대한 일일 종가 자료입니다.
| kos | kq | kl | ki | WonDol | sam | |
|---|---|---|---|---|---|---|
| Date | ||||||
| 2023-01-10 | 2351.310059 | 696.049988 | 14440.0 | 4885.0 | 1239.280029 | 60400.0 |
| 2023-01-11 | 2359.530029 | 709.770020 | 14525.0 | 4875.0 | 1240.160034 | 60500.0 |
| 2023-01-12 | 2365.100098 | 710.820007 | 14580.0 | 4860.0 | 1241.599976 | 60500.0 |
| ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
자료의 각 변수에 대한 VIF를 계산하고 그 결과를 기준으로 선택한 변수들의 VIF를 계산합니다.
다음 코드는 위 자료를 호출하고 정리하기 위한 것입니다.
st=pd.Timestamp(2023,1, 10)
et=pd.Timestamp(2024, 5, 30)
code=["KS11", "KQ11", "122630", "114800","USD/KRW","005930"]
nme=["kos","kq","kl", "ki", "WonDol","sam" ]
da=pd.DataFrame()
for i, j in zip(nme,code):
d=fdr.DataReader(j,st, et)[["Close"]]
da=pd.concat([da, d], axis=1)
da=da.ffill()
da.columns=nme
자료를 표준화합니다.
Scaler=StandardScaler().fit(da) daN=Scaler.transform(da) print(daN[:1, ])
[[-1.87429503 -3.01558934 -1.48676546 1.5873922 -2.36245434 -1.73039517]]
다음의 변수들의 상관계수 결과는 변수 "kos", "kl", "ki" 사이에 매우 밀접한 관계임을 나타냅니다.
data=pd.DataFrame(daN, columns=nme) data.corr().round(3)
| kos | kq | kl | ki | WonDol | sam | |
|---|---|---|---|---|---|---|
| kos | 1.000 | 0.711 | 0.976 | -0.977 | 0.332 | 0.790 |
| kq | 0.711 | 1.000 | 0.592 | -0.598 | 0.302 | 0.480 |
| kl | 0.976 | 0.592 | 1.000 | -0.997 | 0.372 | 0.871 |
| ki | -0.977 | -0.598 | -0.997 | 1.000 | -0.365 | -0.873 |
| WonDol | 0.332 | 0.302 | 0.372 | -0.365 | 1.000 | 0.490 |
| sam | 0.790 | 0.480 | 0.871 | -0.873 | 0.490 | 1.000 |
위 결과는 VIF 결과와 같은 유사합니다.
from statsmodels.stats.outliers_influence import variance_inflation_factor
n, p=daN.shape vif=[variance_inflation_factor(daN, i) for i in range(p)] np.around(pd.DataFrame(vif, index=nme), 4)
| 0 | |
|---|---|
| kos | 153.6790 |
| kq | 7.0381 |
| kl | 222.1241 |
| ki | 236.2386 |
| WonDol | 1.4154 |
| sam | 13.1149 |
위 결과에서"kos","kl","ki"의 VIF는 매우 높습니다. 이 변수들로 다중공선성 문제를 유발할 수 있습니다. 그러므로 세 변수중에서 하나만을 선택할 수 있습니다.
select=daN[:, [1, 3, 4, 5]] nme1=np.array(nme)[[1, 3, 4, 5]] n1, p1=select.shape vif=[variance_inflation_factor(select, i) for i in range(p1)] np.around(pd.DataFrame(vif, index=nme1), 4)
| 0 | |
|---|---|
| kq | 1.6198 |
| ki | 5.3283 |
| WonDol | 1.3837 |
| sam | 5.0493 |
sklearn.linear_model.LinearRegression()의 객체.score() 메서드를 적용한 VIF 계산은 다음과 같습니다.
vif=np.array([])
for i in range(select.shape[1]):
xnew=np.delete(select, i, axis=1)
ynew=select[:,i].reshape(-1,1)
m=LinearRegression().fit(xnew, ynew)
r0=m.score(xnew, ynew)
vif0=1/(1-r0)
vif=np.append(vif, vif0)
re=pd.DataFrame(vif, index=nme1, columns=["VIF"])
re.round(4)
| VIF | |
|---|---|
| kq | 1.6198 |
| ki | 5.3283 |
| WonDol | 1.3837 |
| sam | 5.0493 |
예 3)
예 2의 자료에서 변수 sam를 반응변수로 하여 모든 설명변수를 사용하는 full model과 위 예제에서 vif를 기준으로 선정된 설명변수를 기반으로 하는 모델을 생성합니다. 두 모델을 평가합니다.
자료에서 설명변수와 반응변수를 다음과 같이 분리 합니다.
ind=da.iloc[:-1,:-1] de=da.iloc[1:, -1:] ind.head(1).round(1)
| kos | kq | kl | ki | WonDol | |
|---|---|---|---|---|---|
| Date | |||||
| 2023-01-10 | 2351.3 | 696.0 | 14440.0 | 4885.0 | 1239.39 |
de.head(1)
| sam | |
|---|---|
| Date | |
| 2023-01-11 | 60500.0 |
설명변수를 표준화합니다.
scaler=StandardScaler().fit(ind) indN=scaler.transform(ind) print(indN[0,:].round(3))
[-1.872 -3.012 -1.484 1.585 -2.362]
모델을 평가하기 위해서 모델 생성에 참여하지 않은 데이터를 사용하는 것이 객관적 평가가 될 것입니다. 모델 생성에 참여한 데이터를 훈련데이터(train data)와 참여하지 않은 데이터를 검증데이터(test data)라고 하며 이들 사이의 차이가 작을수록 모델의 일반화에 좋다고 할 수 있습니다. 모델이 검증데이터 보다 훈련데이터에 더 적합한 경우를 과적합(overfitting)이라 하며 반대의 경우를 과소적합(underfitting)이라 합니다. 즉, 모델의 일반화는 이러한 과적합이나 과소적합의 경향이 최소가 되도록 하는 과정입니다.
먼저 훈련과 검증데이터를 생성하기 위해 sklearn.model_selection. train_test_split(x, y, test_size) 함수를 사용합니다.
from sklearn.model_selection import train_test_split
full model을 작성하기 위해 모든 설명변수와 반응변수를 훈련과 검증 데이터로 분리합니다. 생성한 모델을 사용하여 훈련데이터와 검증데이터에 각각 결정계수(R2)를 비교합니다.
Xtr, Xte, ytr, yte=train_test_split(indN, de, test_size=0.2, random_state=3)
m=LinearRegression().fit(Xtr, ytr)
R2_tr=m.score(Xtr, ytr)
R2_te=m.score(Xte, yte)
print("R2_train: %.3f, R2_test: %.3f" %(R2_tr, R2_te))
R2_train: 0.896, R2_test: 0.917
위 결과에 의하면 과소적합이 관찰됩니다. 다음으로 예 2에서 선택한 일부 설명변수에 대 위 과정 동일한 방법으로 R2를 산출합니다.
ind2=ind[["kq", "ki","WonDol"]] scaler2=StandardScaler().fit(ind2) indN2=scaler2.transform(ind2) print(indN2[0])
[-3.01187177 1.58489292 -2.36163584]
Xtr2, Xte2, ytr2, yte2=train_test_split(indN2, de, test_size=0.2, random_state=3)
m2=LinearRegression().fit(Xtr2, ytr2)
R2_tr2=m2.score(Xtr2, ytr2)
R2_te2=m2.score(Xte2, yte2)
print("R_train: %.3f, R_test: %.3f" %(R2_tr2, R2_te2))
R_train: 0.774, R_test: 0.790
위 결과는 full model에 비해 훈련과 검증 데이터 사이의 R2의 차이가 감소됨을 알 수 있습니다. 그러나 결정계수의 수준이 저하됨을 알 수 있습니다. 이것은 변수 선택의 문제가 있음을 의미합니다.
댓글
댓글 쓰기