重回帰分析にダミー変数を使うとき

# 重回帰分析にダミー変数を使うとき

## 1.　重回帰分析に使用するデータ


```python
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` 引数で最初のダミー変数を除くこともできる。


```python
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



```python
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


実際に使用する変数は後で決めることにして，以下のようなデータフレームを作る。


```python
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$ を予測する重回帰式を求める。


```python
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



```python
y = df6['y']
```

以前に作っておいた関数を使って結果を求める。


```python
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))
```


```python
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 個のダミー変数を使う場合


```python
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 個のダミー変数を使う場合


```python
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 個のダミー変数を使う場合


```python
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 個のダミー変数を使った重回帰分析の結果から，どれか一つのダミー変数を除いた結果を導くことができる。


```python
(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$個のダミー変数を使う），予測結果は変わらない。


```python
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]



```python
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() 関数では自動的に一次従属な説明変数はモデルから外される。