重回帰分析にダミー変数を使うとき
1. 重回帰分析に使用するデータ
import pandas as pd
df5 = pd.DataFrame({
'x1': [2, 2, 3, 1, 2, 3, 2, 2, 2, 3, 2, 2, 1, 3, 2, 2, 1, 2, 1, 1],
'x2': [2, 2, 3, 2, 2, 1, 2, 3, 2, 3, 2, 1, 3, 2, 2, 1, 1, 3, 2, 1],
'y': [17, 20, 24, 15, 20, 19, 20, 22, 20, 26, 24, 14, 22, 23, 23, 19, 13, 27, 17, 13]})
前の記事では $x_1$, $x_2$ ともに連続変数として，そのまま説明変数として使った。
今回，$x_2$ をカテゴリー変数として扱うことにする。
pd.get_dummies() により，$k$ 個のカテゴリーを表すために 冗長ではあるが $k$ 個のダミー変数を作る。
$k$ 個のダミー変数は一次従属なので,正規方程式の逆行列を求めて...というような昔ながらのアルゴリズムを使うと，Singular matrix errorが発生し，解を求めることができない。
ちゃんとしたアルゴリズムであれば，$k$ 個のダミー変数を使うようにしても，解は求まる。
しかし，$k$ 個のダミー変数を使った場合と $k-1$ 個のダミー変数を使った場合でどのような違いがあるのか，そもそも違いなんてないのか？
実際に見てみようじゃないかというのが今回のテーマ。
重回帰分析でカテゴリー変数を使用する場合には，どれか 1 つのダミー変数を削除して，$k-1$ 個のダミー変数を使用すればよい。あるいは pd.get_dummies() の drop_first=True 引数で最初のダミー変数を除くこともできる。
dummy3 = pd.get_dummies(df5["x2"], prefix="d")
print(dummy3.head(3))
d_1 d_2 d_3
0 0 1 0
1 0 1 0
2 0 0 1
dummy2 = pd.get_dummies(df5["x2"], prefix="d", drop_first=True)
print(dummy2.head(3))
d_2 d_3
0 1 0
1 1 0
2 0 1
実際に使用する変数は後で決めることにして，以下のようなデータフレームを作る。
df6 = pd.concat([df5, dummy3], axis=1)
print(df6.head(5))
x1 x2 y d_1 d_2 d_3
0 2 2 17 0 1 0
1 2 2 20 0 1 0
2 3 3 24 0 0 1
3 1 2 15 0 1 0
4 2 2 20 0 1 0
2. sklern.linear_model.LinearRegression による分析
2.1. 3 個のダミー変数を使う場合
連続変数 $x1$，ダミー変数 $d_1$，$d_2$，$d_3$ を説明変数として，$y$ を予測する重回帰式を求める。
from sklearn.linear_model import LinearRegression
lr = LinearRegression()
X = df6.drop(['x2', 'y'], axis=1)
print(X.head(3))
x1 d_1 d_2 d_3
0 2 0 1 0
1 2 0 1 0
2 3 0 0 1
y = df6['y']
以前に作っておいた関数を使って結果を求める。
from sklearn.linear_model import LinearRegression
import pandas as pd
def lm(X, y):
lr = LinearRegression()
lr = lr.fit(X, y)
b0, b = lr.intercept_, lr.coef_
beta = tuple(b * X.std(axis=0) / y.std())
R2 = lr.score(X, y)
yhat = lr.predict(X)
Se = sum((y - yhat)**2)
St = sum((y - y.mean())**2)
R2adj = 1 - (Se / (n - p - 1)) / (St / (n - 1))
return(dict(b0=b0, b=b, beta=beta, R2=R2, yhat=yhat, R2adj=R2adj))
X = df6.drop(['x2', 'y'], axis=1)
y = df6['y']
lm(X, y)
{'b0': 14.509019607843138,
'b': array([ 2.74117647, -3.84313725, 0.1827451 , 3.66039216]),
'beta': (0.4571465989589744,
-0.4148666354410395,
0.022779160706274896,
0.3951393035925156),
'R2': 0.773055971922641,
'yhat': array([20.17411765, 20.17411765, 26.39294118, 17.43294118, 20.17411765,
18.88941176, 20.17411765, 23.65176471, 20.17411765, 26.39294118,
20.17411765, 16.14823529, 20.91058824, 22.91529412, 20.17411765,
16.14823529, 13.40705882, 23.65176471, 17.43294118, 13.40705882]),
'R2adj': 0.7125375644353453}
予測式は以下のようになる。
$\hat{y} = b_0 + b_1\ x_1 + b_2\ d_1 + b_3\ d_2 + b_4\ d_3 \tag{1}$
$\hat{y} = 14.509019607843138 + 2.74117647\ x_1 -3.84313725\ d_1 + 0.1827451\ d_2 + 3.66039216\ d_3 \tag{2}$
2.2. 2 個のダミー変数を使う場合
X = df6.drop(['x2', 'y', 'd3_1'], axis=1)
y = df6['y']
results = lm(X, y)
results
{'b0': 10.665882352941175,
'b': array([2.74117647, 4.02588235, 7.50352941]),
'beta': (0.45714659895897425, 0.501825888005632, 0.8100059390335557),
'R2': 0.7730559719226411,
'yhat': array([20.17411765, 20.17411765, 26.39294118, 17.43294118, 20.17411765,
18.88941176, 20.17411765, 23.65176471, 20.17411765, 26.39294118,
20.17411765, 16.14823529, 20.91058824, 22.91529412, 20.17411765,
16.14823529, 13.40705882, 23.65176471, 17.43294118, 13.40705882]),
'R2adj': 0.7125375644353453}
予測値 yhat を見れば判るように，ダミー変数が 3 個の場合も 2 個の場合も同じ予測値になっている。
次式に示すように，$k$ 個のカテゴリーに対する $k$ 個のダミー変数に対する係数は，一次従属なので定数項とのやり取りでどのようにでもなる。例えば，(1) 式と (3) 式は同じになる。
逆にいえば，回帰式を一意的にするには $k-1$ 個のダミー変数を使うようにすれば良い。最も，どのダミー変数を除くかで一意的にはなりえないのではあるが。
$\hat{y} = b_0 + b_1\ x_1 + 0\ d_1 + b_3\ d_2 + b_4\ d_3 \tag{3}$
$\hat{y} = 10.665882352941175 + 2.74117647\ x_1 + 0\ d_1 + 4.02588235\ d_2 + 7.50352941\ d_3 \tag{4}$
$\hat{y} = (10.665882352941175 + 3.84313725) + 2.74117647\ x_1 + (0 - 3.84313725)\ d_1 + (4.02588235 - 3.84313725)\ d_2 + (7.50352941 - 3.84313725)\ d_3 \tag{5}$
$\hat{y} = 14.509019602941175 + 2.74117647\ x_1 - 3.84313725\ d_1 + 0.1827451\ d_2 + 3.66039216\ d_3 \tag{6}$
3. statsmodels.formula.api を用いて重回帰分析を行う
3.1. 3 個のダミー変数を使う場合
import pandas as pd
import statsmodels.formula.api as smf
model = smf.ols(formula='y ~ x1 + d_1 + d_2 + d_3', data=df6)
fitted3 = model.fit();
print(fitted3.summary())
OLS Regression Results
==============================================================================
Dep. Variable: y R-squared: 0.773
Model: OLS Adj. R-squared: 0.731
Method: Least Squares F-statistic: 18.17
Date: Wed, 18 May 2022 Prob (F-statistic): 2.10e-05
Time: 16:44:53 Log-Likelihood: -41.330
No. Observations: 20 AIC: 90.66
Df Residuals: 16 BIC: 94.64
Df Model: 3
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
Intercept 10.8818 1.145 9.504 0.000 8.455 13.309
x1 2.7412 0.733 3.741 0.002 1.188 4.295
d_1 -0.2159 0.810 -0.267 0.793 -1.933 1.501
d_2 3.8100 0.684 5.570 0.000 2.360 5.260
d_3 7.2876 0.939 7.763 0.000 5.297 9.278
==============================================================================
Omnibus: 1.082 Durbin-Watson: 1.850
Prob(Omnibus): 0.582 Jarque-Bera (JB): 0.941
Skew: 0.471 Prob(JB): 0.625
Kurtosis: 2.510 Cond. No. 1.53e+16
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The smallest eigenvalue is 4.72e-31. This might indicate that there are
strong multicollinearity problems or that the design matrix is singular.
$d_1$ の $p$ 値が 0.793 となっており， $d_1$ は使う必要がないことを示唆している。
3.2. 2 個のダミー変数を使う場合
model = smf.ols(formula='y ~ x1 + d_2 + d_3', data=df6)
fitted2 = model.fit();
print(fitted2.summary())
OLS Regression Results
==============================================================================
Dep. Variable: y R-squared: 0.773
Model: OLS Adj. R-squared: 0.731
Method: Least Squares F-statistic: 18.17
Date: Wed, 18 May 2022 Prob (F-statistic): 2.10e-05
Time: 16:45:50 Log-Likelihood: -41.330
No. Observations: 20 AIC: 90.66
Df Residuals: 16 BIC: 94.64
Df Model: 3
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
Intercept 10.6659 1.629 6.549 0.000 7.213 14.119
x1 2.7412 0.733 3.741 0.002 1.188 4.295
d_2 4.0259 1.172 3.434 0.003 1.540 6.511
d_3 7.5035 1.383 5.427 0.000 4.572 10.435
==============================================================================
Omnibus: 1.082 Durbin-Watson: 1.850
Prob(Omnibus): 0.582 Jarque-Bera (JB): 0.941
Skew: 0.471 Prob(JB): 0.625
Kurtosis: 2.510 Cond. No. 9.39
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
この場合も前節と同じく，3 個のダミー変数を使った重回帰分析の結果から，どれか一つのダミー変数を除いた結果を導くことができる。
(10.8818 - 0.2159, -0.2159 + 0.2159, 3.8100 + 0.2159, 7.2876 + 0.2159)
(10.6659, 0.0, 4.0259, 7.503500000000001)
当たり前であるが，$k$個のダミー変数を使おうが,どのダミー変数を除こうが（$k−1$個のダミー変数を使う），予測結果は変わらない。
print(fitted3.predict())
[20.17411765 20.17411765 26.39294118 17.43294118 20.17411765 18.88941176
20.17411765 23.65176471 20.17411765 26.39294118 20.17411765 16.14823529
20.91058824 22.91529412 20.17411765 16.14823529 13.40705882 23.65176471
17.43294118 13.40705882]
print(fitted2.predict())
[20.17411765 20.17411765 26.39294118 17.43294118 20.17411765 18.88941176
20.17411765 23.65176471 20.17411765 26.39294118 20.17411765 16.14823529
20.91058824 22.91529412 20.17411765 16.14823529 13.40705882 23.65176471
17.43294118 13.40705882]
4. 結論
説明変数にカテゴリー変数を使う場合には，「カテゴリー数」個のダミー変数を使うよりは，「カテゴリー数 - 1」個のダミー変数を使うほうがよい。
R の lm() 関数では自動的に一次従属な説明変数はモデルから外される。