분산 분석 기반의 카테고리 분석

회귀 분석 대상이 되는 독립 변수가 카테고리 값을 가지는 변수인 경우에는 카테고리 값에 의해 연속 변수인 y값이 달라진다. 이러한 경우, 분산 분석(ANOVA)을 사용하면 카테고리 값의 영향을 정량적으로 분석할 수 있다. 또한 이는 카테고리 값에 의해 회귀 모형이 달라지는 것으로도 볼 수 있기 때문에 모형 비교에도 사용될 수 있다.

카테고리 독립 변수와 더미 변수

카테고리 값은 여러개의 다른 상태를 나타내는 값이다. 분석시에는 편의상 이 값은 0, 1과 같은 정수로 표현하지만 원래 카테고리값은 1, 2, 3과 같이 숫자로 표현되어 있어도 이는 단지 "A", "B", "C"라는 라벨을 숫자로 대신 쓴 것에 지나지 않으며 실제로 크기의 의미가 없다는 점에 주의해야 한다. 즉, 2라는 값이 1보다 2배 더 크다는 뜻이 아니고 3이라는 값도 1보다 3배 더 크다는 뜻이 아니다.

따라서 카테고리 값을 그냥 정수로 쓰면 회귀 분석 모형은 이 값을 크기를 가진 숫자로 인식할 수 있는 위험이 있기 때문에 반드시 one-hot-encoding 등을 통해 더미 변수(dummy variable)의 형태로 변환해야 함

더미 변수는 0 또는 1만으로 표현되는 값으로 어떤 요인이 존재하는가 존재하지 않는가를 표시하는 독립 변수이다. 다음과 같은 명칭으로도 불린다.

  • indicator variable
  • design variable
  • Boolean indicator
  • binary variable
  • treatment

In [1]:
from sklearn.preprocessing import OneHotEncoder
encoder = OneHotEncoder()

In [2]:
x0 = np.random.choice(3, 10)
x0


Out[2]:
array([2, 2, 1, 0, 2, 0, 2, 0, 0, 0])

In [3]:
encoder.fit(x0[:, np.newaxis])
X = encoder.transform(x0[:, np.newaxis]).toarray()
X


Out[3]:
array([[ 0.,  0.,  1.],
       [ 0.,  0.,  1.],
       [ 0.,  1.,  0.],
       [ 1.,  0.,  0.],
       [ 0.,  0.,  1.],
       [ 1.,  0.,  0.],
       [ 0.,  0.,  1.],
       [ 1.,  0.,  0.],
       [ 1.,  0.,  0.],
       [ 1.,  0.,  0.]])

In [4]:
dfX = pd.DataFrame(X, columns=encoder.active_features_)
dfX


Out[4]:
0 1 2
0 0.0 0.0 1.0
1 0.0 0.0 1.0
2 0.0 1.0 0.0
3 1.0 0.0 0.0
4 0.0 0.0 1.0
5 1.0 0.0 0.0
6 0.0 0.0 1.0
7 1.0 0.0 0.0
8 1.0 0.0 0.0
9 1.0 0.0 0.0

더미 변수와 모형 비교

더미 변수를 사용하면 사실상 회귀 모형 복수개를 동시에 사용하는 것과 실질적으로 동일하다.

더미 변수의 예 1

$$ Y = \alpha_{1} + \alpha_{2} D_2 + \alpha_{3} D_3 $$
  • $D_2 = 0, D_3 = 0$ 이면 $Y = \alpha_{1} $
  • $D_2 = 1, D_3 = 0$ 이면 $Y = \alpha_{1} + \alpha_{2} $
  • $D_2 = 0, D_3 = 1$ 이면 $Y = \alpha_{1} + \alpha_{3} $

더미 변수의 예 2

$$ Y = \alpha_{1} + \alpha_{2} D_2 + \alpha_{3} D_3 + \alpha_{4} X $$
  • $D_2 = 0, D_3 = 0$ 이면 $Y = \alpha_{1} + \alpha_{4} X $
  • $D_2 = 1, D_3 = 0$ 이면 $Y = \alpha_{1} + \alpha_{2} + \alpha_{4} X $
  • $D_2 = 0, D_3 = 1$ 이면 $Y = \alpha_{1} + \alpha_{3} + \alpha_{4} X $

더미 변수의 예 3

$$ Y = \alpha_{1} + \alpha_{2} D_2 + \alpha_{3} D_3 + \alpha_{4} X + \alpha_{5} D_4 X + \alpha_{6} D_5 X $$
  • $D_2 = 0, D_3 = 0$ 이면 $Y = \alpha_{1} + \alpha_{4} X $
  • $D_2 = 1, D_3 = 0$ 이면 $Y = \alpha_{1} + \alpha_{2} + (\alpha_{4} + \alpha_{5}) X $
  • $D_2 = 0, D_3 = 1$ 이면 $Y = \alpha_{1} + \alpha_{3} + (\alpha_{4} + \alpha_{6}) X $

더미 변수의 예 4: Boston Dataset


In [5]:
from sklearn.datasets import load_boston
boston = load_boston()
dfX0_boston = pd.DataFrame(boston.data, columns=boston.feature_names)
dfy_boston = pd.DataFrame(boston.target, columns=["MEDV"])

import statsmodels.api as sm
dfX_boston = sm.add_constant(dfX0_boston)

df_boston = pd.concat([dfX_boston, dfy_boston], axis=1)
df_boston.tail()


Out[5]:
const CRIM ZN INDUS CHAS NOX RM AGE DIS RAD TAX PTRATIO B LSTAT MEDV
501 1 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 1 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 1 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 1 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 1 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

In [6]:
dfX_boston.CHAS.plot()
dfX_boston.CHAS.unique()


Out[6]:
array([ 0.,  1.])

In [9]:
model = sm.OLS(dfy_boston, dfX_boston)
result = model.fit()
print(result.summary())


                            OLS Regression Results                            
==============================================================================
Dep. Variable:                   MEDV   R-squared:                       0.741
Model:                            OLS   Adj. R-squared:                  0.734
Method:                 Least Squares   F-statistic:                     108.1
Date:                Fri, 03 Jun 2016   Prob (F-statistic):          6.95e-135
Time:                        12:42:11   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.4911      5.104      7.149      0.000      26.462      46.520
CRIM          -0.1072      0.033     -3.276      0.001      -0.171      -0.043
ZN             0.0464      0.014      3.380      0.001       0.019       0.073
INDUS          0.0209      0.061      0.339      0.735      -0.100       0.142
CHAS           2.6886      0.862      3.120      0.002       0.996       4.381
NOX          -17.7958      3.821     -4.658      0.000     -25.302     -10.289
RM             3.8048      0.418      9.102      0.000       2.983       4.626
AGE            0.0008      0.013      0.057      0.955      -0.025       0.027
DIS           -1.4758      0.199     -7.398      0.000      -1.868      -1.084
RAD            0.3057      0.066      4.608      0.000       0.175       0.436
TAX           -0.0123      0.004     -3.278      0.001      -0.020      -0.005
PTRATIO       -0.9535      0.131     -7.287      0.000      -1.211      -0.696
B              0.0094      0.003      3.500      0.001       0.004       0.015
LSTAT         -0.5255      0.051    -10.366      0.000      -0.625      -0.426
==============================================================================
Omnibus:                      178.029   Durbin-Watson:                   1.078
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              782.015
Skew:                           1.521   Prob(JB):                    1.54e-170
Kurtosis:                       8.276   Cond. No.                     1.51e+04
==============================================================================

Warnings:
[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.

In [8]:
params1 = result.params.drop("CHAS")
params1


Out[8]:
const      36.491103
CRIM       -0.107171
ZN          0.046395
INDUS       0.020860
NOX       -17.795759
RM          3.804752
AGE         0.000751
DIS        -1.475759
RAD         0.305655
TAX        -0.012329
PTRATIO    -0.953464
B           0.009393
LSTAT      -0.525467
dtype: float64

In [10]:
params2 = params1.copy()
params2["const"] += result.params["CHAS"]
params2


Out[10]:
const      39.179665
CRIM       -0.107171
ZN          0.046395
INDUS       0.020860
NOX       -17.795759
RM          3.804752
AGE         0.000751
DIS        -1.475759
RAD         0.305655
TAX        -0.012329
PTRATIO    -0.953464
B           0.009393
LSTAT      -0.525467
dtype: float64

In [12]:
df_boston.boxplot("MEDV", "CHAS")
plt.show()



In [25]:
sns.stripplot(x="CHAS", y="MEDV", data=df_boston, jitter=True, alpha=.3)
sns.pointplot(x="CHAS", y="MEDV", data=df_boston, dodge=True, color='r')
plt.show()


분산 분석을 이용한 모형 비교

$K$개의 복수의 카테고리 값을 가지는 더미 변수의 영향을 보기 위해서는 F-검정을 통해 복수 개의 모형을 비교하는 분산 분석을 사용할 수 있다.

이 경우에는 분산 분석에 사용되는 각 분산의 의미가 다음과 같다.

  • ESS: 각 그룹 평균의 분산 (Between-Group Variance) $$ BSS = \sum_{k=1}^K (\bar{x} - \bar{x}_k)^2 $$

  • RSS: 각 그룹 내의 오차의 분산의 합 (Within-Group Variance) $$ WSS = \sum_{k=1}^K \sum_{i} (x_{i} - \bar{x}_k)^2 $$

  • TSS : 전체 오차의 분산 $$ TSS = \sum_{i} (x_{i} - \bar{x})^2 $$

source degree of freedom mean square F statstics
Between $$\text{BSS}$$ $$K-1$$ $$\dfrac{\text{ESS}}{K-1}$$ $$F$$
Within $$\text{WSS}$$ $$N-K$$ $$\dfrac{\text{RSS}}{N-K}$$
Total $$\text{TSS}$$ $$N-1$$ $$\dfrac{\text{TSS}}{N-1}$$
$R^2$ $$\text{BSS} / \text{TSS}$$

이 때 F-검정의 귀무가설은 $\text{BSS}=0$ 즉 $\text{WSS}=\text{TSS}$ 이다. 즉, 그룹간 차이가 없는 경우이다.


In [27]:
import statsmodels.api as sm
model = sm.OLS.from_formula("MEDV ~ C(CHAS)", data=df_boston)
result = model.fit()
table = sm.stats.anova_lm(result)
table


Out[27]:
df sum_sq mean_sq F PR(>F)
C(CHAS) 1.0 1312.079271 1312.079271 15.971512 0.000074
Residual 504.0 41404.216144 82.151223 NaN NaN

In [32]:
model1 = sm.OLS.from_formula("MEDV ~ CRIM + ZN +INDUS + NOX + RM + AGE + DIS + RAD + TAX + PTRATIO + B + LSTAT", data=df_boston)
model2 = sm.OLS.from_formula("MEDV ~ CRIM + ZN +INDUS + NOX + RM + AGE + DIS + RAD + TAX + PTRATIO + B + LSTAT + C(CHAS)", data=df_boston)
result1 = model1.fit()
result2 = model2.fit()
table = sm.stats.anova_lm(result1, result2)
table


Out[32]:
df_resid ssr df_diff ss_diff F Pr(>F)
0 493.0 11299.555411 0.0 NaN NaN NaN
1 492.0 11080.276284 1.0 219.279126 9.7367 0.001912