회귀분석 with python

0. Linear Regression

수치를 예측하는 방법 중 하나로 통계학에서는 회귀분석으로 잘 알려져 있다. 학습이 가능하기 때문에 머신러닝으로도 또한 잘 알려져 있다.

단순 선형 회귀와 같은 경우에는 독립변수 하나를 가지고 학습을 하며 이후에 새로운 독립변수들에 대해서 종속변수를 예측하게 된다. 단순 선형 회귀식은 일차 함수처럼 나타나게 된다.

$ \hat{y} = \beta_{1} x + \beta_{0} $

위의 식을 보면 단순한 일차 방정식처럼 보인다는 것을 알 수 있다. 여기서 우리가 구하려는 것은 이미 알고 있는 $x_{i}$ 들과 $y_{i}$들을 토대로 $\beta_{1}$과 $\beta_{0}$를 구하는 것이다.

쉽게 말해 기울기와 y 절편을 구한다고 생각하면 된다.

Linear Regression은 수치 예측에 있어서 너무너무 유명하므로 자세한 설명은 생략하고 코드로 구현하는 것만 작성하였다.
하지만 Multi Linear Regression과 같이 독립변수의 개수가 늘어나게 되고, 또한 이 방법론에서 쓰이는 최소제곱법의 경우에는 기하학적인 이해가 필요하다.
회귀분석의 핵심인 최소제곱법의 기하학적 이해는 따로 포스팅한다 !

1. 데이터 불러오기

데이터는 사이킷런 라이브러리에 있는 보스턴 집값 데이터를 활용하였다. 데이터는 딕셔너리 형태로 받아지므로 key에서 data와 target에 해당하는 부분을 가져와서 데이터 프레임으로 바꾸어 준다.

# 필요한 라이브러리 불러오기
from sklearn.datasets import load_boston
import pandas as pd

house = load_boston()
house.keys()

dict_keys(['data', 'target', 'feature_names', 'DESCR', 'filename'])

- 보스턴 데이터 변수

house['feature_names']

array(['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD',
       'TAX', 'PTRATIO', 'B', 'LSTAT'], dtype='<U7')

- 보스턴 데이터에 대한 설명

house['DESCR'].split("\n")

['.. _boston_dataset:',
 '',
 'Boston house prices dataset',
 '---------------------------',
 '',
 '**Data Set Characteristics:**  ',
 '',
 '    :Number of Instances: 506 ',
 '',
 '    :Number of Attributes: 13 numeric/categorical predictive. Median Value (attribute 14) is usually the target.',
 '',
 '    :Attribute Information (in order):',
 '        - CRIM     per capita crime rate by town',
 '        - ZN       proportion of residential land zoned for lots over 25,000 sq.ft.',
 '        - INDUS    proportion of non-retail business acres per town',
 '        - CHAS     Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)',
 '        - NOX      nitric oxides concentration (parts per 10 million)',
 '        - RM       average number of rooms per dwelling',
 '        - AGE      proportion of owner-occupied units built prior to 1940',
 '        - DIS      weighted distances to five Boston employment centres',
 '        - RAD      index of accessibility to radial highways',
 '        - TAX      full-value property-tax rate per $10,000',
 '        - PTRATIO  pupil-teacher ratio by town',
 '        - B        1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town',
 '        - LSTAT    % lower status of the population',
 "        - MEDV     Median value of owner-occupied homes in $1000's",
 '',
 '    :Missing Attribute Values: None',
 '',
 '    :Creator: Harrison, D. and Rubinfeld, D.L.',
 '',
 'This is a copy of UCI ML housing dataset.',
 'https://archive.ics.uci.edu/ml/machine-learning-databases/housing/',
 '',
 '',
 'This dataset was taken from the StatLib library which is maintained at Carnegie Mellon University.',
 '',
 "The Boston house-price data of Harrison, D. and Rubinfeld, D.L. 'Hedonic",
 "prices and the demand for clean air', J. Environ. Economics & Management,",
 "vol.5, 81-102, 1978.   Used in Belsley, Kuh & Welsch, 'Regression diagnostics",
 "...', Wiley, 1980.   N.B. Various transformations are used in the table on",
 'pages 244-261 of the latter.',
 '',
 'The Boston house-price data has been used in many machine learning papers that address regression',
 'problems.   ',
 '     ',
 '.. topic:: References',
 '',
 "   - Belsley, Kuh & Welsch, 'Regression diagnostics: Identifying Influential Data and Sources of Collinearity', Wiley, 1980. 244-261.",
 '   - Quinlan,R. (1993). Combining Instance-Based and Model-Based Learning. In Proceedings on the Tenth International Conference of Machine Learning, 236-243, University of Massachusetts, Amherst. Morgan Kaufmann.',
 '']

데이터의 독립변수

X = pd.DataFrame(house['data'], columns=house['feature_names'])
X

	CRIM	ZN	INDUS	CHAS	NOX	RM	AGE	DIS	RAD	TAX	PTRATIO	B	LSTAT
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
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
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
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
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
...	...	...	...	...	...	...	...	...	...	...	...	...	...
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
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
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
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
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

506 rows × 13 columns

Y = pd.DataFrame(house['target'], columns=['target'])
Y

	target
0	24.0
1	21.6
2	34.7
3	33.4
4	36.2
...	...
501	22.4
502	20.6
503	23.9
504	22.0
505	11.9

506 rows × 1 columns

 

data와 target 변수를 합친 데이터 프레임을 생성한다.

house = pd.concat([X,Y], axis=1 )

house

	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

2. 독립변수와의 상관관계

pd.DataFrame(house.corr())

	CRIM	ZN	INDUS	CHAS	NOX	RM	AGE	DIS	RAD	TAX	PTRATIO	B	LSTAT	target
CRIM	1.000000	-0.200469	0.406583	-0.055892	0.420972	-0.219247	0.352734	-0.379670	0.625505	0.582764	0.289946	-0.385064	0.455621	-0.388305
ZN	-0.200469	1.000000	-0.533828	-0.042697	-0.516604	0.311991	-0.569537	0.664408	-0.311948	-0.314563	-0.391679	0.175520	-0.412995	0.360445
INDUS	0.406583	-0.533828	1.000000	0.062938	0.763651	-0.391676	0.644779	-0.708027	0.595129	0.720760	0.383248	-0.356977	0.603800	-0.483725
CHAS	-0.055892	-0.042697	0.062938	1.000000	0.091203	0.091251	0.086518	-0.099176	-0.007368	-0.035587	-0.121515	0.048788	-0.053929	0.175260
NOX	0.420972	-0.516604	0.763651	0.091203	1.000000	-0.302188	0.731470	-0.769230	0.611441	0.668023	0.188933	-0.380051	0.590879	-0.427321
RM	-0.219247	0.311991	-0.391676	0.091251	-0.302188	1.000000	-0.240265	0.205246	-0.209847	-0.292048	-0.355501	0.128069	-0.613808	0.695360
AGE	0.352734	-0.569537	0.644779	0.086518	0.731470	-0.240265	1.000000	-0.747881	0.456022	0.506456	0.261515	-0.273534	0.602339	-0.376955
DIS	-0.379670	0.664408	-0.708027	-0.099176	-0.769230	0.205246	-0.747881	1.000000	-0.494588	-0.534432	-0.232471	0.291512	-0.496996	0.249929
RAD	0.625505	-0.311948	0.595129	-0.007368	0.611441	-0.209847	0.456022	-0.494588	1.000000	0.910228	0.464741	-0.444413	0.488676	-0.381626
TAX	0.582764	-0.314563	0.720760	-0.035587	0.668023	-0.292048	0.506456	-0.534432	0.910228	1.000000	0.460853	-0.441808	0.543993	-0.468536
PTRATIO	0.289946	-0.391679	0.383248	-0.121515	0.188933	-0.355501	0.261515	-0.232471	0.464741	0.460853	1.000000	-0.177383	0.374044	-0.507787
B	-0.385064	0.175520	-0.356977	0.048788	-0.380051	0.128069	-0.273534	0.291512	-0.444413	-0.441808	-0.177383	1.000000	-0.366087	0.333461
LSTAT	0.455621	-0.412995	0.603800	-0.053929	0.590879	-0.613808	0.602339	-0.496996	0.488676	0.543993	0.374044	-0.366087	1.000000	-0.737663
target	-0.388305	0.360445	-0.483725	0.175260	-0.427321	0.695360	-0.376955	0.249929	-0.381626	-0.468536	-0.507787	0.333461	-0.737663	1.000000

타겟변수인 집값과 상관관계가 높은 변수는 RM, LSTAT 인 것을 알 수 있다

상관관계가 높아보이는 두 변수를 토대로 산점도를 그려서 선형적인 모습을 띄는지 확인한다.

import matplotlib.pyplot as plt

RM 변수와의 상관관계는 양의 상관관계를 갖는 모습을 볼 수 있다. 즉 RM변수의 값이 커지면 타겟변수의 값도 커지고, RM 변수의값이 작아지면 타겟변수의 값도 작아진다는 말이다.

plt.scatter(house['RM'], house['target'], color='grey')

 

LSTAT변수와의 상관관계는 음의 상관관계를 갖는 모습을 볼 수 있다. 즉 LSTAT 변수의 값이 증가할 때 타겟변수는 감소하고, LSTAT 변수의 값이 감소할 때 타겟변수는 증가하는 관계라 할 수 있다.

plt.scatter(house['LSTAT'], house['target'], color='grey')

3. 회귀분석

# 필요한 라이브러리 불러오기
from sklearn.linear_model import LinearRegression

독립변수와 타겟변수를 명시하여 준다.

single linear regression이기 때문에 독립변수는 하나만 사용하는데 LSTAT 변수를 사용하였다.

x = house[['LSTAT']]
y = house[['target']]

x

	LSTAT
0	4.98
1	9.14
2	4.03
3	2.94
4	5.33
...	...
501	9.67
502	9.08
503	5.64
504	6.48
505	7.88

506 rows × 1 columns

y

	target
0	24.0
1	21.6
2	34.7
3	33.4
4	36.2
...	...
501	22.4
502	20.6
503	23.9
504	22.0
505	11.9

506 rows × 1 columns

회귀 모델 생성하기

회귀 분석 함수를 사용하기 위해 클래스에서 생성자를 호출하여 인스턴스 변수 reg에 담아준다.

reg = LinearRegression(fit_intercept=True)

reg.fit(x,y)

LinearRegression(copy_X=True, fit_intercept=True, n_jobs=None, normalize=False)

생성된 모델, 회귀식 확인하기

$ \beta_{0} $ 의 값을 확인한다

beta_0 = reg.intercept_
beta_0

array([34.55384088])

$ \beta_{1} $의 값을 확인한다.

beta_1 = reg.coef_
beta_1

array([[-0.95004935]])

$ \hat{y} = 34.55x - 0.95 $ 와 같은 회귀식을 생성한 것이다.

이 회귀식은 최소제곱법으로 만들어진 식인데, 이 최소제곱법이라는 것이 선형대수로 계산을 하는 도중에 만족해야할 것들로 인해서 가정들이 있다.

1. 오차항은 정규분포를 따른다.
2. 선형적인 관계를 보인다.
3. 독립변수들은 서로 독립적이다.

4. 모델 학습을 통해 예측하기

모든 케이스에 대하여 랜덤으로 70%, 30%로 분할하여 준다. 70%의 케이스는 훈련데이터로 모델을 학습하고, 30%인 케이스에 대해서는 라벨값을 지워주고 학습한 모델을 토대로 예측을 실시하게 된다.

house.shape

(506, 14)

보스턴 데이터는 케이스가 506개라는 것을 알 수 있다.
훈련데이터와 검증데이터로 나누어주는 좋은 함수가 있다. 불러와서 실행해보자

from sklearn.model_selection import train_test_split

# 훈련/검증 셋으로 나누기
train_features, test_features, train_labels, test_labels = train_test_split(x,y, test_size=0.7)

아까 만들었던 것처럼 모델을 훈련데이터에 대해 만들어 준다.

lm = LinearRegression(fit_intercept=True)
lm.fit(train_features, train_labels)

LinearRegression(copy_X=True, fit_intercept=True, n_jobs=None, normalize=False)

lm.coef_

array([[-0.87412931]])

lm.intercept_

array([33.07000416])

$\hat{y} = 33.07x - 0.87 $ 과 같은 식을 얻었다.

이제 test 셋에 대하여 예측을 해보자.

pred_y = lm.predict(test_features)

pred_y는 훈련데이터로 생성한 모델을 사용하여 검증데이터의 피쳐에 대해 예측을 하여 생성한 예측값 $\hat{y}$ 값들이다.

5. 검증

생성한 모델로 예측한 값들이 실제 값과 얼마나 차이가 나는지를 통해서 이 모델의 성능을 평가 할 수 있다. 대체적으로 RMSE(Root Mean Squared Error)라고 해서 오차 제곱들의 평균에 루트를 씌운 값으로 성능을 측정한다.(수치를 예측하는 경우.)

쉽게 말하면 오차들의 편차를 구한다는 식으로 생각하면 될 것 같다. 당연히 그 편차가 적을 수록 잘 예측했다고 얘기할 수 있는 것이다.

- 실제값과 예측치 데이터 프레임

rmse_df = pd.DataFrame()
rmse_df['y'] = test_labels['target']
rmse_df['predict'] = pred_y
rmse_df

	y	predict
197	30.3	25.543751
308	22.8	29.101457
367	23.1	21.417860
467	19.1	14.433567
359	22.6	21.994786
...	...	...
57	31.6	29.617193
338	20.6	25.631164
238	23.7	27.510542
468	19.1	17.222040
125	21.4	20.124149

355 rows × 2 columns

잔차를 포함한 데이터 프레임

import numpy as np
error = np.power(list(rmse_df['y'] - rmse_df['predict']), 2)
rmse_df['e'] = error

rmse_df

	y	predict	e
197	30.3	25.543751	22.621907
308	22.8	29.101457	39.708361
367	23.1	21.417860	2.829594
467	19.1	14.433567	21.775595
359	22.6	21.994786	0.366284
...	...	...	...
57	31.6	29.617193	3.931522
338	20.6	25.631164	25.312608
238	23.7	27.510542	14.520228
468	19.1	17.222040	3.526735
125	21.4	20.124149	1.627796

355 rows × 3 columns

RMSE를 구하자.

$ \sqrt{ \frac{1}{355} \times \sum_{i=1}^{355} (y_{i} - \hat{y_{i}})^{2}} $

그대로 코드로 입력을 한다.

np.sqrt((1/355) * sum(rmse_df['e']))

6.551441274628164

RMSE가 6.55 가 나왔는데, 이는 오차들의 편차가 6.55라는 의미로 해석이 가능하며 당연히 0에 가까울수록 예측 정확도가 좋다고 말할 수 있다.

6. 결정계수

모델의 성능을 RMSE로 설명을 할 수 있다면, 모델에 대한 설명력을 나타내주는 결정계수라는 것을 구해본다.

결정계수는 $R^{2}$으로 표기한다.

결정계수를 설명하기에 앞서 세가지 Summation에 대해 알아보자.

1. SST(총 변동) : 표본 평균과 실제값 차이의 제곱 합.$\sum_{i=1}^{n}(y_{i} - \bar{y})^{2}$

2. SSR(sum of squares of Regression:회귀식으로 설명된 변동): 표본평균과 예측값의 차이의 제곱 합.$\sum_{i=1}^{n}(\hat{y}_{i} - \bar{y})$

3. SSE(sum of squares of Error:회귀식으로 설명 안되는 변동): 실제값과 예측값의 차이의 제곱 합.$\sum_{i=1}^{n}(y_{i} - \hat{y}_{i})$

총 변동 = 설명 된 변동 + 설명 안된 변동

SST = SSE + SSR

결정계수 $R^{2}$은 총 변동 중에서 회귀식으로 설명이 된 변동의 비율로써 총 변동에 대해 회기식으로 얼마나 변동을 설명할 수 있는지를 나타내는 지표이다.

이를 통해서 결정계수가 높을 수록 이 모델은 변동에 대해 잘 설명하고 있다고 할 수 있기 때문에 더 믿을만하다고 말할 수 있다.

rmse_df

	y	predict	e
197	30.3	25.543751	22.621907
308	22.8	29.101457	39.708361
367	23.1	21.417860	2.829594
467	19.1	14.433567	21.775595
359	22.6	21.994786	0.366284
...	...	...	...
57	31.6	29.617193	3.931522
338	20.6	25.631164	25.312608
238	23.7	27.510542	14.520228
468	19.1	17.222040	3.526735
125	21.4	20.124149	1.627796

355 rows × 3 columns

SST, SSR, SSE 를 구해서 결정계수를 구해보자 !

mean_of_y = np.mean(rmse_df['y'])

SSE = sum(rmse_df['e'])
SSR = sum([np.power(i - mean_of_y,2) for i in rmse_df['predict']])

SST = SSE + SSR

print('결정계수: ', SSR/SST)

결정계수:  0.4897591589594024

즉 결정계수 $R^{2}$이 약 0.48 로써 이 모델은 48% 의 설명력을 가지고 있다고 할 수 있다.

자세히 말하면 총 변동 중 모델로 설명할 수 있는 설명력이 48% 정도라는 것이다. 남은 52% 는 이 모델로 설명할 수 없는 변동이라는 것이다.

마지막으로 수정 결정계수를 구해보자.

유의하지 않은 변수가 추가되어도 증가하는 결정계수를 보완하기 위함.

$R_{adj}^{2} = 1 - [ \frac{n - 1}{n - (p+1)}] \frac{SSE}{SST} $

print('수정된 결정계수: ',1 - ((355-1)/(355-(1+1))) * (SSE/SST))

수정된 결정계수:  0.4883137174833667

p는 변수의 개수라서 수정된 결정계수는 늘 결정계수보다 같거나 작을 수 밖에 없음.

독립변수의 수가 다른 모델을 비교할 때 사용하게 된다.