안녕하세요~ 꽁냥이에요.
여러 개의 변수를 포함하는 데이터를 이용하여 선형 회귀 모형을 적합하는 상황을 생각해봅시다. 이 경우 어떤 변수 집합을 사용할지에 대한 문제가 발생합니다. 변수 집합을 생각하지 않고 모든 변수를 다 써도 되겠지만 1) 많은 변수를 포함하는 모형은 해석이 복잡해질 수 있으며 3) 과적합(Overfitting)의 문제, 4) 변수간 다중공선성(Multicolinearity) 존재 가능성 증가의 문제가 발생할 수 있습니다. 또한 적합한 모형을 계속 사용하기 위해서는 모형에 포함된 변수의 데이터를 유지관리해야 하는데 5) 변수가 많아질 수록 데이터 유지 및 관리가 어려워질 수 있습니다. 따라서 적절한 변수 집합을 선택하는 것이 중요하게 됩니다.
이번 포스팅에서는 최적 모형 선택 방법 중에서 전 모형 탐색(all possible search) 또는 최적 부분집합 알고리즙(Best subsets algorithm)에 대해서 알아보고 실제 데이터를 이용하여 구현해보려고 합니다.
1. 전 모형 탐색법(최적 부분집합 알고리즘)
변수가 $p$개 있는 경우 이를 이용하여 만들 수 있는 모든 모형의 개수는 일차항만 고려할 경우 $2^p$개가 됩니다. 전 모형 탐색법이란 $2^p$개의 모든 모형에 대해서 점수를 구하고 가장 낮은(또는 높은) 점수를 갖는 모형을 찾는 방법입니다.
따라서 변수의 개수 $p$가 증가 할수록 계산해야 되는 양은 엄청나게 많아지는데요. 보통은 $p\leq 30$일 때 사용한다고 하네요.
그렇다면 모형 점수라는 선택기준을 골라야 하는데요. 보통 아래와 같은 5가지를 이용합니다.
수정된 결정계수
Mallow's $C_p$
AIC(Akaike's information criterion)
BIC(Bayes information criterion)
PRESS(PREdiction Sum of Square)
먼저 데이터 $(\tilde{x}_i,y_i), \:\:i=1,2,\cdots,n$가 있다고 할게요. 여기서 $\tilde{x}_i=(x_{i1}, \cdots, x_{ip})^t$ 입니다. 여기서는 최소제곱법을 이용하여 회귀 모형을 적합한다고 하겠습니다. 이렇게 얻어진 모형을 이용한 $i$번째 적합값을 $\hat{y}_i$이라고 하겠습니다. 그리고 전체변동 $SST=\sum_{i=1}^{n}(y_i-\bar{y})^2$, 잔차변동 $SSE=\sum_{i=1}^{n}(y_i-\hat{y}_i)^2$ 입니다.
1. 수정된 결정계수
수정된 결정계수는 기존의 결정계수가 변수가 늘어날수록 커지는 단점을 보완한 것인데요. 다음과 같이 정의합니다.
$$R_{a}^{2} = 1 - \left(\frac{n-1}{n-p}\right)\frac{SSE}{SST}$$
수정된 결정계수는 전체변동($SST$)에 대한 잔차변동($SSE$)의 비율이 포함된다는 점에서 모형 적합이 잘되었는지를 판단하는 척도로 사용됩니다. 따라서 수정된 결정계수 값이 높을수록 좋은 모형이라고 판단합니다.
2. Mallow's $C_p$
Mallow's $C_p$는 다음과 같이 정의합니다.$$C_p = \frac{SSE}{MSE_{full}}-(n-2p)$$여기서 $MSE_{full}$은 모든 변수를 포함시켜 최소제곱모형을 적합시켰을 때 계산된 평균잔차제곱합입니다.$C_p$의 경우 모형의 편의(bias)가 작을수록 $C_p\sim p$를 만족합니다. 따라서 모형 점수를 Mallow's $C_p$로 정하였다면 그 값이 작고 또한 $C_p$값이 변수의 개수$p$와 가까운 값을 갖는 모형을 선호하게 됩니다.
3. AIC(Akaike's information criterion)와 BIC(Bayes information criterion)
AIC와 BIC 또한 모형 선택 기준 척도로써 많이 활용됩니다.
AIC와 BIC는 다음과 같이 정의합니다.
$$AIC = n\log SSE - n\log n + 2p \\ BIC = n\log SSE - n\log n + (n\log n)p$$
두 개 모두 모형 적합 정도를 나타내는 잔차제곱합과 변수의 개수가 많은 모형은 피하도록 벌점(penalty)항으로 이루어져 있습니다. AIC 값이 작은 모형을 더 선호하며 BIC도 마찬가지로 값이 작은 모형을 선호하게 됩니다.
4. PRESS(PREdiction Sum of Square)
PRESS는 다음과 같이 정의합니다.
$$PRESS = \sum_{i=1}^{n}(y_i-\hat{y}_{i(i)})^2$$
여기서 $\hat{y}_{i(i)}$은 $i$ 번째를 제외한 데이터를 이용한 모형 적합값입니다. PRESS는 $i$번째 실제 반응 변수 값과 추정값의 차이를 이용했다는 점에서 예측 오차를 나타내는 척도라 할 수 있습니다. 따라서 이 값이 작을수록 좋은 예측력을 가진 모형이라고 할 수 있습니다.
PRESS를 계산하려면 $n-1$개의 데이터를 이용하여 $n$번의 모형 적합을 해야 하지만 Hat matrix를 사용하면 모든 데이터를 이용하여 적합한 모형 하나로 계산할 수 있습니다. 먼저 모형 행렬(Model matrix) $X = (x_{ij}), \:\:(i=1, \;.\;.\;. , n, \: j=1, \;.\;.\;.,p)$, Hat matrix $H = X(X^tX)^{-1}X^t$라고 한다면 PRESS는 다음과 같습니다.
$$PRESS = \sum_{i=1}^{n}\left(\frac{y_i-\hat{y}_{i}}{1-h_{ii}}\right)^2$$
여기서 $h_{ii}$는 Hat matrix의 $i$번째 대각 원소입니다.
이제 앞서 알아보았던 5가지 기준을 통하여 최적 모형을 선택하는 과정을 구현해볼까요?
2. 전 모형 탐색법 구현
먼저 이번 포스팅에서 필요한 데이터(+데이터 설명 파일)를 다운받아주세요.
필요한 모듈을 임포트하고 데이터를 불러와주세요.
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams['axes.unicode_minus'] = False ## 마이나스 '-' 표시 제대로 출력
from itertools import combinations
from statsmodels.formula.api import ols
df = pd.read_csv('./surgical_unit_reduced.csv') ## 데이터 불러오기
df.insert(0,'Intercept',[1]*len(df)) ## 상수항 추가
이제 각 변수 집합에 대하여 앞서 살펴본 5가지의 통계량을 계산합니다. 변수의 개수가 4개이므로 총 16(=$2^4$)개의 변수 집합이 존재합니다. 꽁냥이는 각 변수 집합에 대하여 5개의 통계량을 계산할 거예요.
fit = ols('Log_ST~Blood_Clotting_Score+Prognostic_Index+Enzyme_Test+Liver_Test',
data=df).fit() ## 풀 모형 적합
## hat matrix 계산
X = np.array(df[['Intercept',
'Blood_Clotting_Score',
'Prognostic_Index',
'Enzyme_Test',
'Liver_Test']]) ## Model matrix for full model
X_tX = np.matmul(X.transpose(),X)
X_tX_inv = np.linalg.inv(X_tX)
hat_matrix = np.matmul(np.matmul(X,X_tX_inv),X.transpose()) ## hat matrix
diagonals = np.array([hat_matrix[i][i] for i in range(len(df))]) ## hat matrix의 대각원소
response = 'Log_ST'
y = df[response] ## 반응 변수 벡터
variables = ['Blood_Clotting_Score',
'Prognostic_Index',
'Enzyme_Test',
'Liver_Test'] ## 총 변수집합
num_var = len(variables) ## 총 변수 개수
num_data = len(df) ## 데이터 개수
mse_full = fit.mse_resid ## 모든 변수를 포함한 mean square residual
mean_response = np.mean(y) ## y의 평균
cp_list = [] ## Mallow's C
ad_r_list = [] ## 수정된 결정계수
aic_list = [] ## Akaike's information criterion
bic_list = [] ## Bayes information criterion
press_list = [] ## press criterion
num_parameter = [] ## 파라미터 즉 절편을 포함한 회귀계수의 개수
subsets = [] ## 변수의 집합
for p in range(num_var+1):
if p == 0: ## 절편만 있는 모형
var_exp = '1'
exp = response + ' ~ ' + var_exp
subsets.append('None') ## 절편만 있고 변수는 없음
sub_fit = ols(exp,data=df).fit() ## 절편만 있는 모형 적합
sse = np.sum(np.square(sub_fit.resid)) ## square sum of residual
sst = np.sum(np.square(y-mean_response)) ## total sum of square
cp = sse/mse_full - (num_data-2*(p+1)) ## Mallow's C
ad_r = 0 ## 수정된 결정계수 절편만 있다면 수정된 결정계수 값은 0
aic = num_data*np.log(sse) - num_data*np.log(num_data) + 2*(p+1) ## Akaike's information criterion
bic = num_data*np.log(sse) - num_data*np.log(num_data) + np.log(num_data)*(p+1) ## Bayes information criterion
press = np.sum(np.square(np.divide(sub_fit.resid,1-diagonals))) ## press
cp_list.append(cp)
ad_r_list.append(ad_r)
aic_list.append(aic)
bic_list.append(bic)
press_list.append(press)
num_parameter.append(p+1)
else:
selected_var = combinations(variables,p)
for s in selected_var:
var_exp = '+'.join(s)
exp = response + ' ~ ' + var_exp
subsets.append(', '.join(s))
sub_fit = ols(exp,data=df).fit()
sse = np.sum(np.square(sub_fit.resid))
sst = np.sum(np.square(y-np.mean(y)))
cp = sse/mse_full - (num_data-2*(p+1))
ad_r = 1 - ((num_data-1)/(num_data-p-1))*(sse/sst)
aic = num_data*np.log(sse) - num_data*np.log(num_data) + 2*(p+1)
bic = num_data*np.log(sse) - num_data*np.log(num_data) + np.log(num_data)*(p+1)
press = np.sum(np.square(np.divide(sub_fit.resid,1-diagonals)))
cp_list.append(cp)
ad_r_list.append(ad_r)
aic_list.append(aic)
bic_list.append(bic)
press_list.append(press)
num_parameter.append(p+1)
df_res = pd.DataFrame()
df_res['Variables'] = subsets
df_res['Number_of_parameter'] = num_parameter
df_res['Ad_R'] = ad_r_list
df_res['Cp'] = cp_list
df_res['AIC'] = aic_list
df_res['BIC'] = bic_list
df_res['PRESS'] = press_list
line 1~2
$C_p$를 계산하기 위하여 변수 4개가 들어간 모형을 적합합니다. 모형 적합 방법은 여기를 참고하세요.
line 5~9
PRESS를 계산하기 위한 Hat matrix와 Hat matrix의 대각 원소를 구하는 코드입니다.
line 36~74
주어진 변수 개수에 대하여 변수의 부분집합을 구합니다. 각 부분 집합으로 이루어진 모형을 적합하여 5개 통계량을 계산합니다.
위 코드를 실행하고 df_res를 확인해볼게요.
변수 개수에 따른 통계량의 변화를 시각적으로 확인해보겠습니다. 코드 설명은 생략할게요.
fig = plt.figure(figsize=(15,15))
fig.set_facecolor('white')
font_size = 12
columns = ['Ad_R', 'Cp', 'AIC', 'BIC', 'PRESS']
xlabel = 'Number of Parameter'
ylabels = ['Adjusted R','C','AIC','BIC','PRESS']
marker_style = dict(color='red',marker='o',markersize=10)
for i in range(len(ylabels)):
ind = 320+i+1
plt.subplot(ind)
plt.scatter(df_res['Number_of_parameter'],df_res[columns[i]],alpha=0.5)
plt.xlabel(xlabel,fontsize=font_size)
plt.ylabel(ylabels[i],fontsize=font_size)
if ylabels[i] == 'Adjusted R':
temp_df = df_res.groupby('Number_of_parameter',as_index=False)[columns[i]].max()
idx = temp_df[columns[i]].idxmax()
max_val = temp_df[columns[i]].max()
plt.scatter(temp_df['Number_of_parameter'][idx], max_val,s=90,color='red')
plt.plot(temp_df['Number_of_parameter'],temp_df[columns[i]])
else:
temp_df = df_res.groupby('Number_of_parameter',as_index=False)[columns[i]].min()
idx = temp_df[columns[i]].idxmin()
min_val = temp_df[columns[i]].min()
plt.scatter(temp_df['Number_of_parameter'][idx], min_val,s=90,color='red')
plt.plot(temp_df['Number_of_parameter'],temp_df[columns[i]])
plt.xticks(range(1,6))
plt.show()
PRESS는 파라미터가 다섯 개(변수의 개수 4 + 절편항)인 모형을 선호하는 반면 나머지 4개의 통계량에 대해서는 파라미터가 4개인 모형을 선호하네요. 이번에는 구체적으로 어떤 변수가 들어간 모형을 추천하는지 확인해보겠습니다.
PRESS 기준으로는 모든 변수를 포함하는 모형을 선택하고 있으며 나머지 4개에 대해서는 Blood_Clotting_Score, Prognostic_Index, Enzyme_Test를 포함한 모형을 선택하고 있습니다.
꽁냥이는 위의 결과와 변수를 절약하는 관점(Parsimony)을 고려하여 최종적으로 Blood_Clotting_Score, Prognostic_Index, Enzyme_Test를 포함한 모형을 선택하였습니다.
모형을 선택할 때에는 하나만을 기준으로 고르기보다는 여러 기준을 이용하여 종합적으로 판단하고 선택하시는 것을 추천합니다.
'데이터 분석 > 데이터 분석' 카테고리의 다른 글
| [회귀 분석] 7. 다중공선성 확인하기 - 분산 팽창 인자 with Python (10) | 2020.10.05 |
|---|---|
| [회귀 분석] 6. 변수 선택법(Variable Selection) with Python (10) | 2020.10.01 |
| [회귀 분석] 4. 오차의 등분산성 검정(테스트)하기 with Python (4) | 2020.09.22 |
| [회귀 분석] 3. 정규분포에 대한 가정 검정하기 with Python (2) | 2020.09.19 |
| [회귀 분석] 2. 잔차도(Residual plot)와 QQ plot 그리기 (0) | 2020.09.18 |
댓글