RAKSUL TechBlog

RAKSULグループのエンジニアが技術トピックを発信するブログです

sklearnとstatsmodelsのデフォルトの回帰モデルで得られる決定係数の違い

TECHNOLOGY Python

この記事は ノバセル Advent Calendar 10日目です。

ノバセルでデータサイエンティストをしている宮野です。前回は CausalImpactのアンチパターン について記載しました。今回はライト目に、自身の備忘録として書き残します。

はじめに

ノバセルでは、比較的シンプルな分析ロジックを用いた分析で意思決定の支援をできるものは、streamlitを活用したWebアプリケーションを社内展開しています。その中には、単回帰分析で得られた決定係数を活用する部分が1部含まれています。

この実装を進める中で、Pythonのsklearnとstatsmodelsを用いて決定係数を求めた際、デフォルト設定で計算される値に差異があることがわかりました。本記事では、この差異の理由を掘り下げていきます。

決定係数について

決定係数とは、回帰モデルにおいて、モデルのあてはまりの良さの判断基準に使われる指標で、以下のように定義されます。

  • y_i ：実測値
  • \hat{y_i} ：予測値
  • \bar{y} ：実測値の平均値

決定係数は、回帰直線のあてはまりがよければ１に近い値をとり、悪ければ０に近い値になります。

sklearnとstatsmodelsで決定係数を求める

実際にサンプルデータを用いて、sklearnとstatsmodelsのデフォルトで求められる決定係数が違うことを示します。

サンプルデータを用意します。

import pandas as pd
import numpy as np
import statsmodels.api as sm
from sklearn.linear_model import LinearRegression
import matplotlib.pyplot as plt


# サンプルデータの生成
np.random.seed(123)
x = np.random.rand(50) * 10
y = 3 * x + np.random.randn(50) * 2

# データフレームの作成
df = pd.DataFrame({'x': x, 'y': y})

得られたデータの散布図が図 1です。

# 散布図の描画
plt.scatter(df['x'], df['y'], alpha=0.7)

# x軸とy軸を0開始に設定
plt.xlim(0, df['x'].max() + 1)
plt.ylim(0, df['y'].max() + 1)

# ラベルを設定
plt.xlabel('x')
plt.ylabel('y')

# グラフを表示
plt.show()

図 1：サンプルデータの散布図

次に、このサンプルデータを用いてsklearnとstatsmodelsで決定係数を求めます。

sklearn

model_lr = LinearRegression()
model_lr.fit(df[['x']], df['y'])
round(model_lr.score(df[['x']], df['y']), 4)

0.9044

statsmodels

model = sm.OLS(df['y'], df[['x']])
res = model.fit()
round(res.rsquared, 4)

0.9825

どちらも求められた決定係数が１に近い値になっていますが、statsmodels で求めたほうがより１に近い値になっています。

極端な例

続いて、極端な例として、回帰モデルのあてはまりが悪くなるようにサンプルデータを作成し、決定係数の違いを確認します。

# サンプルデータの生成
np.random.seed(123)
x_extreme = np.random.rand(50) * 10  # 0から10までのランダムな値
y_extreme = np.random.rand(50) * 10  # 0から10までのランダムな値（x_extremeとは独立）

# データフレームの作成
df_extreme = pd.DataFrame({'x_extreme': x_extreme, 'y_extreme': y_extreme})

# 散布図の描画
plt.scatter(df_extreme['x_extreme'], df_extreme['y_extreme'], alpha=0.7)

# x軸とy軸を0開始に設定
plt.xlim(0, df_extreme['x_extreme'].max() + 1)
plt.ylim(0, df_extreme['y_extreme'].max() + 1)

# ラベルを設定
plt.xlabel('x_extreme')
plt.ylabel('y_extreme')

# グラフを表示
plt.show()

図 2：回帰モデルのあてはまりが悪くなるようなサンプルデータの散布図

sklearn

model_lr.fit(df_extreme[['x_extreme']], df_extreme['y_extreme'])
round(model_lr.score(df_extreme[['x_extreme']], df_extreme['y_extreme']), 4)

0.0679

statsmodels

model_extreme = sm.OLS(df_extreme['y_extreme'], df_extreme[['x_extreme']])
res_extreme = model_extreme.fit()
round(res_extreme.rsquared, 4)

0.5728

散布図からもわかるように、今回作成したサンプルデータでは回帰直線の当てはまりは悪くなる ≒ 決定係数は０に近くなるはずです。

ですが、sklearnで求めた決定係数は０にかなり近いのに対し、statsmodelsで求めた決定係数は０に近いとは言えず、大きく異なる結果になりました。

差分が起こる理由

ここが本題ですが、差分の理由は「切片を含むモデルか、含まないモデルか」です。切片を含まないとは、回帰直線が原点（0, 0）を通ることを意味します。

sklearnのLinearRegression()では、デフォルトでは切片を含むようになっているのに対し、statsmodels のデフォルトでは切片を含まないようになっています。この違いが決定係数の差分を引き起こしています。

極端な例で用意したサンプルデータの散布図（図 2）に、切片を含む回帰直線と、切片を含まない回帰直線を引いた図が以下の通りです。

図 3：切片を含む回帰直線
図 4：切片を含まない（原点を通る）回帰直線

切片を含まないモデルで得られる決定係数は、回帰モデルのあてはまりを表す指標ではなくなるため注意が必要です。

切片を含まない場合の決定係数は以下のように定義されます。

statsmodelsで切片を含むモデルの決定係数を得るには、以下のように定数項を追加すると、sklearnで得られる決定係数と一致します。

# 定数項を追加
x_extreme_with_const = sm.add_constant(df_extreme[['x_extreme']])

model_intercept_extreme = sm.OLS(df_extreme['y_extreme'], x_extreme_with_const)
res_intercept_extreme = model_intercept_extreme.fit()
round(res_intercept_extreme.rsquared, 4)

0.0679

おわりに

生成AIの普及により、データサイエンティストなどの専門家に聞かなくとも、コードや数式をすぐに生成し結果が得られる時代です。しかし、本記事で示したように、sklearnとstatsmodelsでは、デフォルト設定のわずかな違いが、分析結果に基づく意思決定や施策の方向性に影響を及ぼす可能性があります。

ツールやライブラリを用いる際には、それらのデフォルト設定や実装の違いを正確に理解することが求められます。

参考資料

回帰直線が原点（0, 0）を通る場合の決定係数については「Rで学ぶ確率統計学　多変量統計編」の第3章にまとめられています。