交互作用項／多項式項の導入前の中心化
さて、化学研究者のモデリング第一歩である線形回帰、その中でも非線形性を簡単に導入できるが故に使用場面の多い応答曲面法についてですが、高次項（交互作用項、多項式項）を導入する前に主効果の中心化を行っていますか？
結論から言うと、筆者は中心化を行った方が良いと考えます。
理由は２つ。
- 多重共線性の程度を定量的に評価するための指標であるVIFが大きくなるのを防ぐ
- 切片や係数の値の解釈性が上がる
本記事では上記２つの理由について、トイデータを用いて解説していきます。
モデリングにおける一番の関心である回帰結果（精度や係数値）は、実は中心化有無で全く変わらない、というのもお示ししていきます。
※筆者自身まだまだ勉強中の身ですので、誤った内容を含んでおりましたらご指摘いただけますと幸いです！
線形回帰や応答曲面法に関する詳説は以下の記事を参照ください。
線形回帰はたった2行で実装できる 化学研究者が最も一般的に用いるモデリング手法は何を隠そう「線形回帰」です。 本記事では、この線形回帰をPythonでどのように実装するかを解説します。 実際に手を動かしながら学習できるようにサンプ[…]
化学分野で活用範囲の広い応答曲面モデル 化学分野でよく使用されるモデリングに線形重回帰がありますが、線形重回帰に最も簡単に非線形性を導入する方法として応答曲面モデルが存在します。 今回は応答曲面モデルについて、Pythonでの実装方[…]
実験用トイデータの作成と確認
データを作ろう
データサイエンスを学習するうえで、自分でトイデータを作ってとにかく実験してみる、というのが非常に重要です。
トイデータの作成から実験まで、Pythonで簡単に実施できるので楽しいですね。
さて、トイデータを作っていきましょう。
今回はx1, x2, x3の３つの説明変数だとします。そして、真の関数は以下の通りに設定します。
x3は全く目的変数に影響しないということになります。
最後のεは誤差項で、適当にデータをばらつかせます。
以下がデータ生成のコードです。
データ点数や誤差項の大きさ、乱数シードを設定していろいろなデータを作れるように、これらを引数にした関数にしています。
（※本記事のサンプルコードではimportを都度行います。これはどの作業でどのライブラリが使われているかを明示するためです。実際はプログラムの最初にまとめて行うのがベター！）
import numpy as np
import pandas as pd
# データ生成の関数
def create_data(
n_samples: int,
noise_std: int,
random_seed: int = 100,
):
# 乱数シードを設定
np.random.seed(random_seed)
# 説明変数の生成
x1 = np.random.uniform(0.1, 1.5, n_samples)
x2 = np.random.uniform(100, 1000, n_samples)
x3 = np.random.uniform(20, 150, n_samples)
# 応答変数 y の定義（20〜95%になるよう設計）
y = (
30
+ 100 * x1
- 40 * (x1**2)
+ 0.1 * x2
- 0.1 * x1 * x2
+ np.random.normal(0, noise_std, n_samples)
)
# 数値を丸める
x1_rounded = np.round(x1, 2)
x2_rounded = np.round(x2).astype(int)
x3_rounded = np.round(x3).astype(int)
y_rounded = np.round(y, 1)
# データフレーム作成
df = pd.DataFrame(
{
"x1": x1_rounded,
"x2": x2_rounded,
"x3": x3_rounded,
"y": y_rounded,
}
)
return dfでは、この関数を使ってデータを生成しましょう。
以降、本記事で同じ結果を得るために同じ設定でデータを生成することをお勧めします。
# 関数を使ったデータ生成
df = create_data(30, 5) # 乱数シードはデフォルト
データを確認しよう
作成したデータセットの外観を確認しましょう。
データフレームでデータ生成しているのでdescribe()関数で要約統計量（平均、標準偏差など）を取得してみましょう。
# 作成したデータの確認
df.describe()
無事に30点のデータが関数で設定した範囲で生成されていることが分かります。
次に、散布図行列を描画してデータのばらつきを確認してみましょう。
matplotlibとseabornを使います。
import matplotlib.pyplot as plt
import seaborn as sns
# 散布図行列を描画
sns.pairplot(df, corner=True, plot_kws={"alpha": 0.6})
plt.suptitle("Pairplot", fontsize=16, y=1.02)
plt.tight_layout()
plt.show()
良い感じにばらついていて、x1, x2, x3には多重共線性の心配もなさそうです。
中心化によるVIFへの影響
まずは理由の１つ目であるVIFについて議論してきましょう。
VIFとは
多重共線性（説明変数間に強い相関関係があること）は線形回帰を行うときの厄介な問題です。
多重共線性があると、モデルの解釈が困難になったり、予測精度が低下する可能性があります。
多重共線性については本ブログでも別記事で詳説予定です。
そんな厄介な多重共線性の程度を定量的に評価するための指標がVIFです。
VIFのイメージは「ある説明変数を他のすべての説明変数でどのくらい説明できるか」です。
数学的な定義は割愛しますが、このVIFは1以上の値を取り、一般的には10を超えると強い多重共線性が疑われる、とされます。
以降で、主効果の中心化をせずに高次項（交互作用項、多項式項）を導入すると、VIFが非常に大きくなってしまうことを確かめます。
中心化無しのVIF
中心化無しのVIFが大きくなることは、勘のいい人は当たり前だと感じるかもしれません。
高次項は当然主効果を使って算出されます。そのため、もちろん主効果で説明できてしまいますよね。
実際にVIFを確かめてみましょう。
まずは中心化をせずに高次項を導入していってみましょう。
コードは以下の通りです。
最後に、高次項を導入した後の散布図行列を描画しています。
import itertools
# 目的変数の列名
target_column = "y"
# 説明変数の列名リストを自動取得（目的変数以外）
explanatory_vars = df.drop(columns=target_column).columns.tolist()
# 中心化無しのデータセットの作成
# 回帰用データフレームの作成
df_no_centering = df.copy()
# 二次項の追加
for var in explanatory_vars:
df_no_centering[f"{var}^2"] = df_no_centering[var] ** 2
# 交互作用項の追加
for var1, var2 in itertools.combinations(explanatory_vars, 2):
df_no_centering[f"{var1}*{var2}"] = df_no_centering[var1] * df[var2]
# 散布図行列を描画
sns.pairplot(df_no_centering, corner=True, plot_kws={"alpha": 0.6})
plt.suptitle("Pairplot", fontsize=16, y=1.02)
plt.tight_layout()
plt.show()
ここでは散布図行列の一部を載せましたが、赤枠で囲った部分で説明変数同士の相関が見えます。
とくに２次項と主効果の相関が非常に強く出ています。
さて、次にVIFを算出しましょう。
VIFの算出にはstatsmodelsのモジュールを使います。
from statsmodels.stats.outliers_influence import variance_inflation_factor
# 説明変数の列名リストを自動取得（目的変数以外）
explanatory_vars_2 = df_no_centering.drop(columns=target_column).columns.tolist()
# 説明変数と目的変数に分ける
X = df_no_centering.drop(columns=target_column).to_numpy()
y = df_no_centering[target_column].to_numpy().reshape(-1, 1)
# 各列のVIFを算出
vif_no_centering = pd.DataFrame()
vif_no_centering["Variable"] = explanatory_vars_2
vif_no_centering["VIF"] = [variance_inflation_factor(X, i) for i in range(X.shape[1])]
print(vif_no_centering) # 結果の出力
この通り、いずれの説明変数のVIFも10を超えてしまっております。
これは強い多重共線性が疑われる水準です。
中心化ありのVIF
では、主効果を中心化してから高次項を導入するとどうなるでしょう。
先ほどと同様にやっていきましょう。
# 中心化有りのデータセットの作成
# 回帰用データフレームの作成
df_centering = df.copy()
# 主効果の中心化
for var in explanatory_vars:
df_centering[var] = df_centering[var] - df_centering[var].mean()
# 二次項の追加
for var in explanatory_vars:
df_centering[f"{var}^2"] = df_centering[var] ** 2
# 交互作用項の追加
for var1, var2 in itertools.combinations(explanatory_vars, 2):
df_centering[f"{var1}*{var2}"] = df_centering[var1] * df_centering[var2]
# 散布図行列を描画
sns.pairplot(df_centering, corner=True, plot_kws={"alpha": 0.6})
plt.suptitle("Pairplot", fontsize=16, y=1.02)
plt.tight_layout()
plt.show()
散布図行列ではぱっと見で直線的な相関は無くなりました。
続いてVIFも計算しましょう。
# 説明変数と目的変数に分ける
X_centering = df_centering.drop(columns=target_column).to_numpy()
y_centering = df_centering[target_column].to_numpy().reshape(-1, 1)
# 各列のVIFを算出
vif_centering = pd.DataFrame()
vif_centering["Variable"] = explanatory_vars_2
vif_centering["VIF"] = [
variance_inflation_factor(X_centering, i) for i in range(X_centering.shape[1])
]
print(vif_centering)
VIFがすべて10以下に落ち着きました。
これが中心化の効果の一つ目です。
高いVIFは回帰結果に影響するのか？
さて、中心化をせずに高次項を導入するとVIFが高くなり、多重共線性があるデータセットになってしまうことが確認されました。
多重共線性があると予測精度が低下する可能性があると知られていますが、この場合もそうなのでしょうか。
結論、影響しないです。
後ほど実際に確かめていきますが、回帰結果には影響しません。
少し調べたところでは、これは最小二乗法の解法による影響だと考えています。
最も基本的な最小二乗法の解法は正規方程式の解を求める方法ですが、これは逆行列を直接計算するので多重共線性があると非常に不安定になります。
一方、今回用いたstatsmodels.apiのOLS.fit()や、scikit-learnのLinearRegressionのfitメソッドは行列分解に基づく解法で、計算の不安定性を改善しているようです。
そのため、今回の実験では回帰結果には影響しなかったと考えられます。
つまり、中心化によってVIFを小さくしておくのは本質的には意味がないです。
ただ、見た目でVIFが大きすぎるのは気持ち悪いので、というくらいでしょうか。
では、なぜ筆者は中心化を推奨するのか。
それは次にお話しする係数の解釈性です。
中心化と係数解釈性
ここで取り扱う内容を簡単にまとめます。
まず、ここまでで中心化有無の違いで作成したデータセットを使って、回帰を行います。
回帰結果として、R2値と係数値、p値を確認します。
そして、最後に係数解釈性について考えてみたいと思います。
まずは回帰してみよう
今回はstatsmodels.apiのOLS.fit()を使います。
scikit-learnのLinearRegressionと違い、回帰結果として係数のp値も返してくれるのでこちらを使用します。
以下がpythonコードです。
回帰の前に標準化を実施しています。
import statsmodels.api as sm
from sklearn.preprocessing import StandardScaler
# 中心化無しデータセットのモデリング
# スケーラーの準備とフィット
X_scaler = StandardScaler()
y_scaler = StandardScaler()
X_scaled = X_scaler.fit_transform(X)
y_scaled = y_scaler.fit_transform(y).ravel() # yは1次元に戻す
# 定数項の追加（OLSでは必須）
X_scaled_const = sm.add_constant(X_scaled)
# statsmodelsで回帰
model = sm.OLS(y_scaled, X_scaled_const)
results = model.fit()
# 中心化ありデータセットのモデリング
X_centering_scaler = StandardScaler()
y_centering_scaler = StandardScaler()
X_centering_scaled = X_centering_scaler.fit_transform(X_centering)
y_centering_scaled = y_centering_scaler.fit_transform(y_centering).ravel()
X_centering_scaled_const = sm.add_constant(X_centering_scaled)
model_centering = sm.OLS(y_centering_scaled, X_centering_scaled_const)
results_centering = model_centering.fit()中心化有無のそれぞれの回帰結果をresults / results_centeringという変数に格納しました。
モデル精度の確認
モデル精度としてR2値を比較してみましょう。
statsmodels.apiのOLS.fit()で回帰した場合、R2値は回帰結果の属性値として格納されているので以下コードのように簡単にアクセスできます。
print(f"R-squared(without centering): {results.rsquared:.3f}")
print(f"R-squared(with centering): {results_centering.rsquared:.3f}")いずれも0.919程度で、中心化有無で精度は全く変わっていません。
係数値（標準化スケール）の比較
次に係数値やそれぞれのp値を比較していきましょう。
ここは可視化すると比較しやすいので、matplotlibを使って棒グラフで表現してみましょう。
import matplotlib.pyplot as plt
# 各モデルの標準化スケールでの係数値とp値の可視化
# ラベル設定
labels = explanatory_vars_2.copy()
labels.insert(0, "intercept")
# 横軸の位置
x = np.arange(len(labels)) # ラベルのインデックス
width = 0.25 # 各棒の幅
# グラフを2つ並べる準備
fig, axs = plt.subplots(2, 2, figsize=(12, 8))
# 中心化無しの回帰結果
axs[0, 0].bar(x - width, results.params, width=width)
axs[0, 0].set_title("No Centering")
axs[0, 0].set_ylabel("Value")
axs[0, 0].axhline(0, color="gray", linestyle="--")
axs[0, 0].set_xticks(x)
axs[0, 0].set_xticklabels(labels, rotation=45, ha="right")
for i, v in enumerate(results.params):
axs[0, 0].text(i, v, f"{v:.2f}", ha="center", va="bottom")
axs[0, 1].bar(x - width, results.pvalues, width=width)
axs[0, 1].set_title("No Centering")
axs[0, 1].set_ylabel("p-value")
axs[0, 1].set_xticks(x)
axs[0, 1].set_xticklabels(labels, rotation=45, ha="right")
for i, v in enumerate(results.pvalues):
axs[0, 1].text(i, v, f"{v:.2f}", ha="center", va="bottom")
# 中心化ありの回帰結果
axs[1, 0].bar(x - width, results_centering.params, width=width)
axs[1, 0].set_title("Centering")
axs[1, 0].set_ylabel("Value")
axs[1, 0].axhline(0, color="gray", linestyle="--")
axs[1, 0].set_xticks(x)
axs[1, 0].set_xticklabels(labels, rotation=45, ha="right")
for i, v in enumerate(results_centering.params):
axs[1, 0].text(i, v, f"{v:.2f}", ha="center", va="bottom")
axs[1, 1].bar(x - width, results_centering.pvalues, width=width)
axs[1, 1].set_title("Centering")
axs[1, 1].set_ylabel("p-value")
axs[1, 1].set_xticks(x)
axs[1, 1].set_xticklabels(labels, rotation=45, ha="right")
for i, v in enumerate(results_centering.pvalues):
axs[1, 1].text(i, v, f"{v:.2f}", ha="center", va="bottom")
# 表示
plt.tight_layout()
plt.show()
中心化有無に依らず、いずれもx1, x2, x1^2, x1*x2の標準化スケールでの係数値が大きく、これらのp値が低い結果となっており、真の関数の特徴を捉えていることになります。
係数の解釈性
標準化スケールでの係数値は相対的にどの説明変数が寄与しているのかの比較には使えますが、実際に各説明変数がどの程度目的変数に影響するかは元のスケールの係数値を確認する必要があります。
そこで、まずは各係数を元のスケールに戻してみましょう。
目的変数及び説明変数をそれぞれ標準化しており、標準化は「平均を引いて標準偏差で割る」という定義を考慮すれば、以下コードの数式の意味は簡単に理解できると思います。
分からない人は実際の数式を書いてみてください。
# 中心化無しの係数、切片のスケールを戻す
X_mean = X_scaler.mean_
X_scale = X_scaler.scale_
y_mean = y_scaler.mean_[0]
y_scale = y_scaler.scale_[0]
coef = results.params[1:] * (y_scale / X_scale)
intercept = y_mean - np.sum(coef * X_mean) + results.params[0]
# 切片も係数リストに追加
coef_list = coef.tolist()
coef_list.insert(0, intercept)
print(coef_list) # 結果の出力
# 中心化ありの係数、切片のスケールを戻す
X_centering_mean = X_centering_scaler.mean_
X_centering_scale = X_centering_scaler.scale_
y_centering_mean = y_centering_scaler.mean_[0]
y_centering_scale = y_centering_scaler.scale_[0]
coef_centering = results_centering.params[1:] * (y_centering_scale / X_centering_scale)
intercept_centering = (
y_centering_mean
- np.sum(coef_centering * X_centering_mean)
+ results_centering.params[0]
)
coef_list_centering = coef_centering.tolist()
coef_list_centering.insert(0, intercept_centering)
print(coef_list_centering) # 結果の出力さて、スケールを戻した係数を使って、それぞれ数式に戻してみましょう。
ここでは、p値から優位だと判断できるx1, x2, x1^2, x1*x2だけで表現します。
まずは中心無しのケースです。以下数式になります。
今回は実験なので真の関数を知っています。
そのため、真の関数と比較しやすいと思われるかもしれませんが、実際には真の関数形を知っていることはないです。
では中心化ありの数式を見てみましょう。
当然中心化をしているので真の関数との比較はできません。でも、問題はそこではないですね。
では、解釈性について考えてみましょう。
分かりやすいのは切片の解釈です。
中心化が無い場合の切片はどう解釈できますか？これには意味を見出せません。
一方、中心化をしている場合、切片はx1, x2が平均値の時のyの値になるので意味を持ちます。
係数値についてはどうでしょう。
中心化している場合、平均から離れていったときにどう目的変数に影響するかを意味します。
このように中心化していると係数の解釈性が増すという利点があります。
中心化有無で回帰結果は変わっているのか？
最後に回帰結果が変わっているかの議論に戻りましょう。
精度は変わっていないことをR2値から確認しました。
最後に各説明変数の係数値が変わっていないことを確認しましょう。
上述の通り、中心化をする場合の数式は上記の通りなので、カッコを展開して中心化無しの時の数式と同じ形に変更してみましょう。
コードは以下の通りです。
# 中心化前の係数、切片を計算
mean_x1 = df["x1"].mean()
mean_x2 = df["x2"].mean()
mean_x3 = df["x3"].mean()
# 係数の計算
coef_original = [
coef_list_centering[1] - 2 * coef_list_centering[4] * mean_x1 - coef_list_centering[7] * mean_x2 - coef_list_centering[8] * mean_x3,
coef_list_centering[2] - 2 * coef_list_centering[5] * mean_x2 - coef_list_centering[7] * mean_x1 - coef_list_centering[9] * mean_x3,
coef_list_centering[3] - 2 * coef_list_centering[6] * mean_x3 - coef_list_centering[8] * mean_x1 - coef_list_centering[9] * mean_x2,
coef_list_centering[4],
coef_list_centering[5],
coef_list_centering[6],
coef_list_centering[7],
coef_list_centering[8],
coef_list_centering[9],
]
# 切片の計算
intercept_original = (
intercept_centering
- coef_list_centering[1] * mean_x1
- coef_list_centering[2] * mean_x2
- coef_list_centering[3] * mean_x3
+ coef_list_centering[4] * (mean_x1**2)
+ coef_list_centering[5] * (mean_x2**2)
+ coef_list_centering[6] * (mean_x3**2)
+ coef_list_centering[7] * mean_x1 * mean_x2
+ coef_list_centering[8] * mean_x1 * mean_x3
+ coef_list_centering[9] * mean_x2 * mean_x3
)
# 切片も係数リストに追加
coef_original.insert(0, intercept_original)
print(coef_original) # 結果の出力
print(coef_list) # 中心化無しの回帰結果の結果を再度出力出力される値が全く同じことが確認できたかと思います。
つまり、回帰結果は中心化有無で変わらないということになります。
おわりに
今回、高次項（交互作用項、多項式項）を導入した応答曲面法を用いる場合に、主効果を中心化する効果についてまとめました。
statsmodelsやscikit-learnに実装されている最小二乗法を用いる場合、中心化効果は回帰結果に影響しません。
一方、求まった切片や係数値が解釈しやすいという利点があるため、筆者は中心化をお勧めします。