[통계] 다중회귀분석 예제 - Statsmodel을 이용한 고유값, vif 확인

다중회귀분석 모델링¶

사용할 데이터셋: sklearn의 보스턴 집값 데이터셋

• 1970년도 인구조사로부터 가져온 보스턴의 506개 인구 조사 구역으로 구성
• 21개의 특성 변수 포함
• 목표 변수 - 주택의 중앙값(median)

In [1]:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib as mpl
%matplotlib inline

from sklearn.datasets import load_boston
from sklearn import linear_model

In [2]:

dataset = load_boston()
boston = pd.DataFrame(dataset.data, columns=dataset.feature_names)
boston['target'] = dataset.target

In [3]:

boston

Out[3]:

	CRIM	ZN	INDUS	CHAS	NOX	RM	AGE	DIS	RAD	TAX	PTRATIO	B	LSTAT	target
0	0.00632	18.0	2.31	0.0	0.538	6.575	65.2	4.0900	1.0	296.0	15.3	396.90	4.98	24.0
1	0.02731	0.0	7.07	0.0	0.469	6.421	78.9	4.9671	2.0	242.0	17.8	396.90	9.14	21.6
2	0.02729	0.0	7.07	0.0	0.469	7.185	61.1	4.9671	2.0	242.0	17.8	392.83	4.03	34.7
3	0.03237	0.0	2.18	0.0	0.458	6.998	45.8	6.0622	3.0	222.0	18.7	394.63	2.94	33.4
4	0.06905	0.0	2.18	0.0	0.458	7.147	54.2	6.0622	3.0	222.0	18.7	396.90	5.33	36.2
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
501	0.06263	0.0	11.93	0.0	0.573	6.593	69.1	2.4786	1.0	273.0	21.0	391.99	9.67	22.4
502	0.04527	0.0	11.93	0.0	0.573	6.120	76.7	2.2875	1.0	273.0	21.0	396.90	9.08	20.6
503	0.06076	0.0	11.93	0.0	0.573	6.976	91.0	2.1675	1.0	273.0	21.0	396.90	5.64	23.9
504	0.10959	0.0	11.93	0.0	0.573	6.794	89.3	2.3889	1.0	273.0	21.0	393.45	6.48	22.0
505	0.04741	0.0	11.93	0.0	0.573	6.030	80.8	2.5050	1.0	273.0	21.0	396.90	7.88	11.9

506 rows × 14 columns

In [4]:

# target을 제외한 column을 변수로 두고 X, y 정의
variables = boston.columns[:-1]
X = boston.iloc[:, :-1]
y = boston['target'].values
size = len(boston)
print("Size of boston dataset: {}".format(size))

Size of boston dataset: 506

Statsmodel 패키지로 모델 구축¶

In [5]:

import statsmodels.api as sm
import statsmodels.formula.api as smf

In [6]:

Xc = sm.add_constant(X)
lr = sm.OLS(y, Xc)
model = lr.fit()

In [7]:

model.summary()

Out[7]:

OLS Regression Results

Dep. Variable:	y	R-squared:	0.741
Model:	OLS	Adj. R-squared:	0.734
Method:	Least Squares	F-statistic:	108.1
Date:	Mon, 06 Sep 2021	Prob (F-statistic):	6.72e-135
Time:	10:32:44	Log-Likelihood:	-1498.8
No. Observations:	506	AIC:	3026.
Df Residuals:	492	BIC:	3085.
Df Model:	13
Covariance Type:	nonrobust

	coef	std err	t	P>|t|	[0.025	0.975]
const	36.4595	5.103	7.144	0.000	26.432	46.487
CRIM	-0.1080	0.033	-3.287	0.001	-0.173	-0.043
ZN	0.0464	0.014	3.382	0.001	0.019	0.073
INDUS	0.0206	0.061	0.334	0.738	-0.100	0.141
CHAS	2.6867	0.862	3.118	0.002	0.994	4.380
NOX	-17.7666	3.820	-4.651	0.000	-25.272	-10.262
RM	3.8099	0.418	9.116	0.000	2.989	4.631
AGE	0.0007	0.013	0.052	0.958	-0.025	0.027
DIS	-1.4756	0.199	-7.398	0.000	-1.867	-1.084
RAD	0.3060	0.066	4.613	0.000	0.176	0.436
TAX	-0.0123	0.004	-3.280	0.001	-0.020	-0.005
PTRATIO	-0.9527	0.131	-7.283	0.000	-1.210	-0.696
B	0.0093	0.003	3.467	0.001	0.004	0.015
LSTAT	-0.5248	0.051	-10.347	0.000	-0.624	-0.425

Omnibus:	178.041	Durbin-Watson:	1.078
Prob(Omnibus):	0.000	Jarque-Bera (JB):	783.126
Skew:	1.521	Prob(JB):	8.84e-171
Kurtosis:	8.281	Cond. No.	1.51e+04

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.51e+04. This might indicate that there are
strong multicollinearity or other numerical problems.

1. $ R^2 $값이 0.741로, 1에 가까운 편이니까 독립변수 X들이 target 값을 잘 표현하고 있다고 할 수 있다(일반적으로 Bio는 0.95, 공학은 0.7, 사회과학은 0.3 정도를 기준으로 한다고 함). 근데 Adjusted $ R^2 $은 무엇일까...? 표준 $ R^2 $과 Adjusted $ R^2 $ 측정값 사이의 비율이 20%를 초과하면 모델에 중복 변수가 있음을 뜻한다고 함. Adjusted $ R^2 $을 어떻게 만드는지는 모르겠음..

2. p-value가 0.5를 넘는 변수들은 제거해도 될 것 같다.

3. Notes에서 "The condition number is large, 1.51e+04. This might indicate that there are strong multicollinearity or other numerical problems."라고 경고하고 있다. 조건수가 통상적으로 30을 넘으면 불안정한 결과로 신뢰도가 낮아진다고 한다. 이유는 다중공선성 때문.

In [8]:

# 상관관계 행렬을 통해 다중공선성을 유발하는 변수들 살펴보기
corr_matrix = X.corr()
print(corr_matrix)

             CRIM        ZN     INDUS      CHAS       NOX        RM       AGE  \
CRIM     1.000000 -0.200469  0.406583 -0.055892  0.420972 -0.219247  0.352734   
ZN      -0.200469  1.000000 -0.533828 -0.042697 -0.516604  0.311991 -0.569537   
INDUS    0.406583 -0.533828  1.000000  0.062938  0.763651 -0.391676  0.644779   
CHAS    -0.055892 -0.042697  0.062938  1.000000  0.091203  0.091251  0.086518   
NOX      0.420972 -0.516604  0.763651  0.091203  1.000000 -0.302188  0.731470   
RM      -0.219247  0.311991 -0.391676  0.091251 -0.302188  1.000000 -0.240265   
AGE      0.352734 -0.569537  0.644779  0.086518  0.731470 -0.240265  1.000000   
DIS     -0.379670  0.664408 -0.708027 -0.099176 -0.769230  0.205246 -0.747881   
RAD      0.625505 -0.311948  0.595129 -0.007368  0.611441 -0.209847  0.456022   
TAX      0.582764 -0.314563  0.720760 -0.035587  0.668023 -0.292048  0.506456   
PTRATIO  0.289946 -0.391679  0.383248 -0.121515  0.188933 -0.355501  0.261515   
B       -0.385064  0.175520 -0.356977  0.048788 -0.380051  0.128069 -0.273534   
LSTAT    0.455621 -0.412995  0.603800 -0.053929  0.590879 -0.613808  0.602339   

              DIS       RAD       TAX   PTRATIO         B     LSTAT  
CRIM    -0.379670  0.625505  0.582764  0.289946 -0.385064  0.455621  
ZN       0.664408 -0.311948 -0.314563 -0.391679  0.175520 -0.412995  
INDUS   -0.708027  0.595129  0.720760  0.383248 -0.356977  0.603800  
CHAS    -0.099176 -0.007368 -0.035587 -0.121515  0.048788 -0.053929  
NOX     -0.769230  0.611441  0.668023  0.188933 -0.380051  0.590879  
RM       0.205246 -0.209847 -0.292048 -0.355501  0.128069 -0.613808  
AGE     -0.747881  0.456022  0.506456  0.261515 -0.273534  0.602339  
DIS      1.000000 -0.494588 -0.534432 -0.232471  0.291512 -0.496996  
RAD     -0.494588  1.000000  0.910228  0.464741 -0.444413  0.488676  
TAX     -0.534432  0.910228  1.000000  0.460853 -0.441808  0.543993  
PTRATIO -0.232471  0.464741  0.460853  1.000000 -0.177383  0.374044  
B        0.291512 -0.444413 -0.441808 -0.177383  1.000000 -0.366087  
LSTAT   -0.496996  0.488676  0.543993  0.374044 -0.366087  1.000000  

In [9]:

corr = np.corrcoef(X, rowvar=0)
eigenvalues, eigenvectors = np.linalg.eig(corr)

In [10]:

print(eigenvalues)

[6.12684883 1.43327512 1.24261667 0.85757511 0.83481594 0.65740718
 0.53535609 0.39609731 0.06350926 0.27694333 0.16930298 0.18601437
 0.22023782]

0에 가까운 고유값은 다중공선성의 요인임을 나타낸다. 따라서 가장 작은 값에 대한 고유벡터를 추출해보자. 다른 값들과 현저히 다른 값의 고유벡터를 가진 변수는 제거할 수 있다.

In [11]:

print(eigenvectors[: ,8])

[-0.0459523   0.08091897  0.25107654 -0.03592171 -0.04363045 -0.0455671
  0.03855068  0.01829854  0.63348972 -0.72023345 -0.02339805  0.00446307
 -0.02443168]

In [12]:

print(variables[2], variables[8], variables[9])

INDUS RAD TAX

In [13]:

# 독립변수들의 VIF 체크
from statsmodels.stats.outliers_influence import variance_inflation_factor
vif = pd.DataFrame()
vif["VIF Factor"] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
vif["features"] = X.columns
vif

Out[13]:

	VIF Factor	features
0	2.100373	CRIM
1	2.844013	ZN
2	14.485758	INDUS
3	1.152952	CHAS
4	73.894947	NOX
5	77.948283	RM
6	21.386850	AGE
7	14.699652	DIS
8	15.167725	RAD
9	61.227274	TAX
10	85.029547	PTRATIO
11	20.104943	B
12	11.102025	LSTAT

보통 VIF 10 이상이면 다중공선성이 요인이 된다고 판단한다는데, 이 데이터에서는 VIF 10 이상인 독립변수가 너무... 많다.. 일단 고유벡터값이 명백하게 다른 INDUS RAD TAX는 모두 10이 넘는다. 그 외에 VIF가 85;;인 PTRATIO와 77인 RM을 추가적으로 삭제해보기로 했다.

추가적으로 히트맵도 그려서 판단해봤다(근데 PTRATIO가 다른 변수들과의 상관관계가 0.5를 안넘어서;; 왜 VIF가 높게 나오는지는 모르겠지만 일단 고유벡터 기준으로 3개의 변수를 먼저 삭제하기로 함)

In [14]:

def make_heatmap(data):
    R = np.corrcoef(data, rowvar=0)  # array 형식 - shape(506, 506)
    R[np.where(np.abs(R) < 0.5)] = 0.0
    heatmap = plt.pcolor(R, cmap=mpl.cm.Blues, alpha=0.8)
    heatmap.axes.set_frame_on(False)
    heatmap.axes.set_xticks(np.arange(R.shape[1]) + 0.5)
    heatmap.axes.set_yticks(np.arange(R.shape[0]) + 0.5)
    heatmap.axes.set_xticklabels(variables)
    heatmap.axes.set_yticklabels(variables)
    plt.xticks(rotation=90)
    plt.colorbar()
    plt.show()
    
make_heatmap(X)

변수선택¶

In [15]:

new_data_1 = X.drop(['INDUS', 'RAD', 'TAX'], axis=1)

In [16]:

Dc_1 = sm.add_constant(new_data_1)
lr_1 = sm.OLS(y, Dc_1)
model_1 = lr_1.fit()

In [17]:

model_1.summary()

Out[17]:

OLS Regression Results

Dep. Variable:	y	R-squared:	0.729
Model:	OLS	Adj. R-squared:	0.723
Method:	Least Squares	F-statistic:	133.1
Date:	Mon, 06 Sep 2021	Prob (F-statistic):	2.35e-133
Time:	10:33:18	Log-Likelihood:	-1510.0
No. Observations:	506	AIC:	3042.
Df Residuals:	495	BIC:	3089.
Df Model:	10
Covariance Type:	nonrobust

	coef	std err	t	P>|t|	[0.025	0.975]
const	29.4163	4.885	6.022	0.000	19.818	39.014
CRIM	-0.0615	0.030	-2.022	0.044	-0.121	-0.002
ZN	0.0415	0.014	3.061	0.002	0.015	0.068
CHAS	3.0430	0.870	3.497	0.001	1.334	4.753
NOX	-15.7814	3.371	-4.682	0.000	-22.404	-9.158
RM	4.1761	0.416	10.037	0.000	3.359	4.994
AGE	-0.0044	0.013	-0.325	0.745	-0.031	0.022
DIS	-1.4510	0.198	-7.330	0.000	-1.840	-1.062
PTRATIO	-0.8366	0.118	-7.112	0.000	-1.068	-0.605
B	0.0084	0.003	3.097	0.002	0.003	0.014
LSTAT	-0.5193	0.051	-10.089	0.000	-0.620	-0.418

Omnibus:	192.122	Durbin-Watson:	1.054
Prob(Omnibus):	0.000	Jarque-Bera (JB):	937.157
Skew:	1.617	Prob(JB):	3.16e-204
Kurtosis:	8.830	Cond. No.	9.29e+03

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 9.29e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

In [18]:

# 다시 독립변수들의 VIF 체크
vif = pd.DataFrame()
vif["VIF Factor"] = [variance_inflation_factor(new_data_1.values, i) for i in range(new_data_1.shape[1])]
vif["features"] = new_data_1.columns
vif

Out[18]:

	VIF Factor	features
0	1.727211	CRIM
1	2.666541	ZN
2	1.132455	CHAS
3	56.732022	NOX
4	73.363974	RM
5	21.293993	AGE
6	13.355022	DIS
7	70.491350	PTRATIO
8	18.453971	B
9	10.996396	LSTAT

In [19]:

## 아직도 PTRATIO랑 RM의 VIF값이 너무 크다. 제거해주자.
new_data_2 = new_data_1.drop(['PTRATIO', 'RM'], axis=1)

In [20]:

Dc_2 = sm.add_constant(new_data_2)
lr_2 = sm.OLS(y, Dc_2)
model_2 = lr_2.fit()

In [21]:

model_2.summary()

Out[21]:

OLS Regression Results

Dep. Variable:	y	R-squared:	0.632
Model:	OLS	Adj. R-squared:	0.626
Method:	Least Squares	F-statistic:	106.8
Date:	Mon, 06 Sep 2021	Prob (F-statistic):	7.07e-103
Time:	10:34:52	Log-Likelihood:	-1587.1
No. Observations:	506	AIC:	3192.
Df Residuals:	497	BIC:	3230.
Df Model:	8
Covariance Type:	nonrobust

	coef	std err	t	P>|t|	[0.025	0.975]
const	44.7672	2.981	15.018	0.000	38.910	50.624
CRIM	-0.0946	0.035	-2.723	0.007	-0.163	-0.026
ZN	0.1015	0.015	6.907	0.000	0.073	0.130
CHAS	4.1261	1.005	4.104	0.000	2.151	6.102
NOX	-15.7653	3.864	-4.080	0.000	-23.358	-8.173
AGE	0.0190	0.015	1.243	0.215	-0.011	0.049
DIS	-1.7912	0.228	-7.854	0.000	-2.239	-1.343
B	0.0058	0.003	1.871	0.062	-0.000	0.012
LSTAT	-0.8825	0.049	-18.096	0.000	-0.978	-0.787

Omnibus:	120.784	Durbin-Watson:	1.121
Prob(Omnibus):	0.000	Jarque-Bera (JB):	245.271
Skew:	1.299	Prob(JB):	5.50e-54
Kurtosis:	5.209	Cond. No.	6.90e+03

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 6.9e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

• 결론: 조건수가 많이 낮아진 것을 보아 다중공선성 문제를 조금 완화한 것 같긴 하다. 근데 변수를 무작정 삭제하는 것만이 방법은 아닐텐데.. $ R^2 $이 같이 줄어든다(그래봤자 -1.0 정도지만)