R/Python の線形回帰では完全な多重共線性のある入力データにどう対処しているか？

この記事は？

R/Python の代表的な統計・機械学習パッケージの線形回帰において、回帰係数推定に用いられている最適化手法と入力データの説明変数間に完全な多重共線性がある場合に見られる挙動を比較した。

サマリー

線形回帰において回帰係数の最小二乗推定量は、計画行列 $\boldsymbol{X}$ と目的変数ベクトル $\boldsymbol{Y}$ を用いて $\hat{\boldsymbol{\beta}} = (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{Y}$ のような形でよく書き表される。

しかし、説明変数間に完全な多重共線性がある場合は $\boldsymbol{X}^\top \boldsymbol{X}$ が逆行列を持たないため上記の式を適用できない。

一方、統計ソフトウェアの内部では、回帰係数の推定の際に上記の式をそのまま用いず、最適化手法として QR分解や擬逆行列を用いたものが採用されていることが多い。

以下の表は、R/Python の代表的なパッケージ3種類について、用いられている最適手法と入力データに完全な多重共線性がある場合に見せる挙動をまとめたものである：

言語	パッケージ	最適化手法	完全な多重共線性のあるデータに対する挙動
R	stats (lm 関数)	QR分解	説明変数を減らして推定を行う（summary に注意書きあり）
Python	statsmodels (OLS)	擬逆行列(default)	最小ノルム推定値を返す（summary に注意書きあり
		QR分解	不安定な推定値を返す （summary に注意書きあり）
Python	scikit-learn (LinearRegression)	擬逆行列	最小ノルム推定値を返す （warningや注意書きなし）

前提

設定

$n$ 組の目的変数 $y_i$ と $p-1$ 個の説明変数 $x_{1i}, x_{2i}, ..., x_{p-1,i}$ からなるデータセット($i=1, 2, ..., n$)をもとに、線形回帰モデルを作成する：

$$y_i = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + ... \beta_{p-1} x_{p-1,i} + \varepsilon_i.$$

ただし、$\varepsilon_i$ は独立同時な分布に従う誤差項とし、その期待値は $E[\varepsilon_i] = 0$、分散は $V[\varepsilon] = \sigma^2$ とする。

この時上記のモデル式は、目的変数ベクトル $\boldsymbol{Y} = (y_1, y_2, ..., y_n)^\top$、計画行列 $\boldsymbol{X}$、回帰係数ベクトル $\boldsymbol{\beta} = (\beta_0, \beta_1, \beta_2, ..., \beta_{p-1})^\top$、誤差項ベクトル $\boldsymbol{\varepsilon} = (\varepsilon_1, \varepsilon_2, ..., \varepsilon_n)^\top$ を用いて以下のように表すことができる：

$$\boldsymbol{Y} = \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$$

ただし、計画行列 $\boldsymbol{X} \in \mathbb{R}^{n \times p}$ は以下のように定義される：

$$\boldsymbol{X} = \begin{pmatrix} 1 & x_{11} & x_{21} & ... & x_{p-1,1}\\ 1 & x_{12} & x_{22} & ... & x_{p-1,2}\\ \vdots & \vdots & \vdots & ... & \vdots\\ 1 & x_{1n} & x_{2n} & ... & x_{p-1,n}\\ \end{pmatrix}$$

最小二乗法による係数推定

最小二乗法では、以下に示す残差二乗和が最小になるように回帰係数を推定する：

$$|| \boldsymbol{Y} - \boldsymbol{X} \boldsymbol{\beta}||^2 = \sum_{i=1}^n (y_i - \beta_0 - \beta_1 x_{1i} - \beta_2 x_{2i} - ... - \beta_{p-1} x_{p-1,i})^2$$

それにより得られる最小二乗推定量 $\hat{\boldsymbol{\beta}}$ は、標準的な教科書^[1]では以下のように紹介されることが多い：

$$\hat{\boldsymbol{\beta}} = \left( \boldsymbol{X}^\top \boldsymbol{X} \right)^{-1} \boldsymbol{X}^\top \boldsymbol{Y} \tag{1}$$

ただし上記の式は、計画行列 $\boldsymbol{X}$ が階数 $p$ でフルランク (${\rm rank}\, \boldsymbol{X} = p$) の時のみ成立する。

もし ${\rm rank}\, \boldsymbol{X} < p$ である場合、$\boldsymbol{X}^\top \boldsymbol{X} \in \mathbb{R}^{p \times p}$ が逆行列を持たないため、上記の式の形では最小二乗推定量を表すことができない（推定量が存在しない or 求まらないわけではない）。
そのような事象が見られるケースとしては、サンプルサイズ $n$ がパラメータ数 $p$ より少ない場合と、説明変数間に完全な多重共線性がある場合が挙げられる。

完全な多重共線性とは？

完全な多重共線性とは、ある説明変数が常に別の説明変数の線型結合で表せることを指す。
例えば、説明変数 $x_{1i}$ と $x_{2i}$ が比例関係（例： $x_{1i} = 3 x_{2i}$）にあったり、または $x_{4i} = 1 - x_{1i} - x_{2i} - x_{3i}$ のような関係が成り立つ場合をいう（後者の例は、カテゴリカル変数を One-Hot Encoding した際によく見られる）。

このようなケースでは、計画行列 $\boldsymbol{X}$ の列ベクトル同士が線型独立ではなくなるため、$\boldsymbol{X}$ はフルランクでなくなる（${\rm rank}\, \boldsymbol{X} < p$）。

多重共線性については、下記の記事も参照いただきたい：

https://zenn.dev/tatamiya/articles/eeb97aae490f3ee8f0e9

最小二乗法における最適化手法

線形回帰係数の最小二乗推定量は計画行列 $\boldsymbol{X}$ が full-rank な場合は式(1)で表すことができたが、実際に計算機で推定量を求める際は大抵の場合は直接この式を適用しない。

以下では、統計ソフトウェアでよく使われる最適化手法として、QR 分解と擬逆行列(pseudo-inverse)を用いた2通りの方法を紹介する。

QR 分解による方法

行列 $\boldsymbol{X} \in \mathbb{R}^{n \times p}$ （ただし、$p\le n$）は、以下のように直交行列 $\boldsymbol{Q} \in \mathbb{R}^{n \times p}$ （$\boldsymbol{Q}^\top \boldsymbol{Q} = \boldsymbol{I}_p$）と上三角行列 $\boldsymbol{R} \in \mathbb{R}^{p \times p}$ の積で表すことができる：

$$\boldsymbol{X} = \boldsymbol{Q}\boldsymbol{R}.$$

これを用いると、

$$\begin{aligned} \boldsymbol{X}^\top \boldsymbol{X} &= \boldsymbol{R}^\top \boldsymbol{Q}^\top \boldsymbol{Q} \boldsymbol{R}\\ &= \boldsymbol{R}^\top \boldsymbol{R} \end{aligned}$$

であることから、回帰係数推定量の式(1)は以下のように書き表すことができる：

$$\begin{aligned} \hat{\boldsymbol{\beta}} &= \left( \boldsymbol{R}^\top \boldsymbol{R} \right)^{-1} \boldsymbol{R}^\top \boldsymbol{Q}^\top \boldsymbol{Y}\\\ &= \boldsymbol{R}^{-1} \boldsymbol{Q}^\top \boldsymbol{Y} \end{aligned}$$

なお、この方法では $\boldsymbol{X}$ が full-rank でない時は逆行列 $\boldsymbol{R}^{-1}$ を求めることができず上記の形で推定量を表すことができない。

擬逆行列を用いた方法

行列 $\boldsymbol{A}$ が与えられたとき、以下のような性質を持つ行列 $\boldsymbol{A}^+$ を擬逆行列 (pseudo-inverse) もしくは Moore-Penrose の逆行列と呼ぶ：

1. $\boldsymbol{A}\boldsymbol{A}^+ \boldsymbol{A} = \boldsymbol{A}$
2. $\boldsymbol{A}^+ \boldsymbol{A} \boldsymbol{A}^+ = \boldsymbol{A}^+$
3. $(\boldsymbol{A} \boldsymbol{A}^+)^\top = \boldsymbol{A} \boldsymbol{A}^+$
4. $(\boldsymbol{A}^+ \boldsymbol{A})^\top = \boldsymbol{A}^+ \boldsymbol{A}$

https://zenn.dev/tatamiya/articles/13df806ad8af9aa457ce

計画行列 $\boldsymbol{X}$ の 擬逆行列 $\boldsymbol{X}^+$ を用いることで、回帰係数 $\boldsymbol{\beta}$ の最小二乗推定量 $\hat{\boldsymbol{\beta}}$ は、以下のように表すことができる：

$$\hat{\boldsymbol{\beta}} = \boldsymbol{X}^+ \boldsymbol{Y}.$$

https://zenn.dev/tatamiya/articles/ca2ab1f8d3f069b78a82

上記は計画行列 $\boldsymbol{X}$ が full-rank の時は式(1)と一致するほか、full-rank でない場合も無数に取りうる推定量の中でノルムが最小になるものと一致する。

パッケージごとの比較

ここでは、以下の3種類の異なるパッケージを対象に、最小二乗法により線形回帰モデルを作成する際にどのような最適化手法が用いられ、完全な多重共線性のあるデータの入力に対してどのような結果を返すかを調べた：

• R: デフォルトで入っている stats パッケージの lm 関数
• Python: statsmodels の OLS クラス
• Python: scikit-learn の LinearRegression クラス

なお、サンプルとして以下のように $x_2 = 3 x_1$ という多重共線性を持つデータセット $\{(x_{1i}, x_{2i}, y_i)\}_{i=1}^5$ を用いた：

（R の場合）

x1 <- c(30, 20, 15, 45, 25)
x2 <- x1 * 3
y <- 2 + 0.3 * x1 + 0.1 * rnorm(5)

df <- data.frame(x1, x2, y)

（Python の場合）

import pandas as pd
import numpy as np

# データ準備
df = pd.DataFrame({"x1": [30, 20, 15, 45, 25]})
df["x2"] = df["x1"] * 3
df["y"] = 2 + 0.3 * df["x1"] +  0.1 * np.random.randn(df.shape[0])

R: stats パッケージの lm 関数の場合

https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/lm

R の lm 関数では、QR 分解を用いて最小二乗推定量を計算する。

ただし、入力説明変数に完全な多重共線性があった場合は、変数を減らして fit した結果を返した。
またこのとき、サマリーに (1 not defined because of singularities) という注意書きが表示された。

res <- lm(y~x1+x2, data=df)
summary(res)

Call:
lm(formula = y ~ x1 + x2, data = df)

Residuals:
       1        2        3        4        5
 0.00140 -0.09743  0.07571  0.01315  0.00718

Coefficients: (1 not defined because of singularities)
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.062194   0.090078   22.89 0.000183 ***
x1          0.296699   0.003117   95.18 2.56e-06 ***
x2                NA         NA      NA       NA
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.07177 on 3 degrees of freedom
Multiple R-squared:  0.9997,	Adjusted R-squared:  0.9996
F-statistic:  9059 on 1 and 3 DF,  p-value: 2.557e-06

なお、lm 関数には singular.ok というオプション引数がある（デフォルトでは TRUE）。
これを FALSE にすると、以下のように入力データに完全な多重共線性がある場合はエラーとなり実行を停止する：

res <- lm(y~x1+x2, data=df, singular.ok = FALSE)

Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :
  singular fit encountered

補足： 完全な多重共線性があった場合どのように変数選択を行うか？

上記で見たように、lm 関数では入力説明変数に完全な多重共線性が認められた場合は変数を減らして推定を行なっていた。

この際に除外する変数は、QR 分解を行う過程で決定される。

QR 分解を実行する際には、入力行列（この場合は計画行列 $\boldsymbol{X}$ ）の先頭から1列ずつ処理を適用していく。
このとき、ある変数列のノルムの大きさが一定値^[2]以下になった場合、その変数は除外対象となる。^[3]

なおそのため、入力データの説明変数の並び順を入れ替えると、除外される変数も変わる場合がある。

Python: statsmodels の OLS の場合

https://www.statsmodels.org/dev/generated/statsmodels.regression.linear_model.OLS.html

statsmodels の OLS では、fit 時の method オプションで最適化手法を以下の2つから選択できる：

• "pinv" (default) ... 擬逆行列によるノルム最小推定量を返す
• "qr" ... QR 分解を用いて求めた推定量を返す

いずれの場合でも、入力説明変数に完全な多重共線性があればサマリー末尾の Notes: に以下のような注意書きが表示された：

The smallest eigenvalue is 2.02e-28. This might indicate that there are
strong multicollinearity problems or that the design matrix is singular.

なお、R の lm 関数の場合と異なり、QR 分解を用いた際に変数選択は行わない。
そのため、係数推定値が極端に大きな（もしくは小さな）値となり、さらに標準誤差がもとまらないなど、不安定さが窺える結果が得られた。

method = "pinv" (default)

import statsmodels.api as sm

# モデル作成（推定方法： Moore-Penrose の逆行列）
mod = sm.OLS(df["y"], sm.add_constant(df[["x1", "x2"]]))
res = mod.fit()

print(res.summary())

        OLS Regression Results
==============================================================================
Dep. Variable:                      y   R-squared:                       1.000
Model:                            OLS   Adj. R-squared:                  1.000
Method:                 Least Squares   F-statistic:                 1.824e+04
Date:                Tue, 13 Feb 2024   Prob (F-statistic):           8.95e-07
Time:                        14:27:43   Log-Likelihood:                 8.9634
No. Observations:                   5   AIC:                            -13.93
Df Residuals:                       3   BIC:                            -14.71
Df Model:                           1
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          2.0057      0.065     30.721      0.000       1.798       2.213
x1             0.0305      0.000    135.062      0.000       0.030       0.031
x2             0.0915      0.001    135.062      0.000       0.089       0.094
==============================================================================
Omnibus:                          nan   Durbin-Watson:                   2.383
Prob(Omnibus):                    nan   Jarque-Bera (JB):                0.365
Skew:                          -0.555   Prob(JB):                        0.833
Kurtosis:                       2.278   Cond. No.                     1.44e+16
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The smallest eigenvalue is 2.02e-28. This might indicate that there are
strong multicollinearity problems or that the design matrix is singular.
/usr/local/lib/python3.10/dist-packages/statsmodels/stats/stattools.py:74: ValueWarning: omni_normtest is not valid with less than 8 observations; 5 samples were given.
  warn("omni_normtest is not valid with less than 8 observations; %i "

method = "qr"

回帰係数推定値は $(-2.622 \times 10^{12}, 8.741 \times 10^{11})$ と極端な値を取り、さらに標準誤差も nan が返ってきた。

res_qr = mod.fit(method="qr")
print(res_qr.summary())

                            OLS Regression Results
==============================================================================
Dep. Variable:                      y   R-squared:                       1.000
Model:                            OLS   Adj. R-squared:                  1.000
Method:                 Least Squares   F-statistic:                 1.400e+04
Date:                Tue, 13 Feb 2024   Prob (F-statistic):           1.33e-06
Time:                        14:26:42   Log-Likelihood:                 8.3018
No. Observations:                   5   AIC:                            -12.60
Df Residuals:                       3   BIC:                            -13.38
Df Model:                           1
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          1.9985      0.075     26.816      0.000       1.761       2.236
x1         -2.622e+12        nan        nan        nan         nan         nan
x2          8.741e+11        nan        nan        nan         nan         nan
==============================================================================
Omnibus:                          nan   Durbin-Watson:                   1.447
Prob(Omnibus):                    nan   Jarque-Bera (JB):                0.287
Skew:                          -0.536   Prob(JB):                        0.866
Kurtosis:                       2.522   Cond. No.                     1.44e+16
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The smallest eigenvalue is 2.02e-28. This might indicate that there are
strong multicollinearity problems or that the design matrix is singular.
/usr/local/lib/python3.10/dist-packages/statsmodels/stats/stattools.py:74: ValueWarning: omni_normtest is not valid with less than 8 observations; 5 samples were given.
  warn("omni_normtest is not valid with less than 8 observations; %i "
/usr/local/lib/python3.10/dist-packages/statsmodels/regression/linear_model.py:1884: RuntimeWarning: invalid value encountered in sqrt
  return np.sqrt(np.diag(self.cov_params()))

Python: scikit-learn の LinearRegression の場合

https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html

scikit-learn の LinearRegression では、最小ノルム推定量（擬逆行列を用いて求めたものと同様のもの）が返ってくる。^[4]
なお、注意書き・警告等の表示は一切行われなかった。

from sklearn.linear_model import LinearRegression

reg = LinearRegression()
reg.fit(df[["x1", "x2"]], df["y"])

print(reg.intercept_)
print(reg.coef_)

2.005712241700081
[0.03051598 0.09154795]

なお、この推定値は statsmodels でデフォルトオプションの擬逆行列 (method="pinv")を用いた場合の値と一致していた。

まとめ

以上のように、R/Python の代表的な統計・機械学習パッケージの線形回帰では、回帰係数の推定に QR 分解か擬逆行列のいずれかが用いられていた。

そして、入力説明変数間に完全な多重共線性がある場合の応答は、アルゴリズムによって異なっていた。
警告・注意書き等が一切表示されない場合や、不安定な推定値が返ってくる場合もあるため、注意が必要である。

なお、実際の分析の過程で完全な多重共線性が見つかった場合の対応であるが、背景となる原因が明確な場合は手動で除外する変数を決めるのが無難である。^[5]
擬逆行列を用いた方法でも確かに推定量は求まるものの、解釈のしやすさの点で扱いが難しくなる。

脚注

1. 例えば、佐和隆光 回帰分析(新装版) (統計ライブラリー、朝倉書店、2020) など

2. 着目する列ベクトルの初期状態におけるノルムの大きさの $10^{-7}$ 倍

3. 以下に調査メモを残している。機会があればこちらも記事の形に整理して公開したい。https://zenn.dev/tatamiya/scraps/7031df482d47fe

4. 詳述すると、直接的には SciPy の scipy.linalg.lstsq 関数を呼び出している。さらにその内部を辿ると、LAPACK (Liner Algegra PACKage) の dgelsd 関数に行きつく。次のメモも参照： https://zenn.dev/link/comments/e1cf0dfab50888

5. あくまでも完全な多重共線性の場合の話である。完全でない多重共線性が認められた場合は、安易に変数を除外するとかえって推定結果にバイアスが生じることがある。