化学分野で活用範囲の広い応答曲面モデル
化学分野でよく使用されるモデリングに線形重回帰がありますが、線形重回帰に最も簡単に非線形性を導入する方法として応答曲面モデルが存在します。
今回は応答曲面モデルについて、Pythonでの実装方法から、作成したモデルの解釈方法まで解説していきます。
化学実験を想定した架空のデータセットを使って解説していくので、ご自身の研究を思い浮かべながら読み進めていくと、具体的に活用方法がイメージしやすいかもしれません！
線形重回帰の基本は以下の記事を参考にしてください。
「線形回帰」とは？～線形単回帰と線形重回帰～ 実験計画法を用いたデータ駆動型の研究開発は、企業研究者が最短で目標達成するのに必要なスキルです。 データ駆動型の研究プロセスの中で、花形とも言える「モデリング」には様々な手法があります。 […]
また、モデルに非線形性を導入する別の手法については以下の記事を参考にしてください。
非線形な関係をモデリングするには 化学分野で使用頻度の高いモデルである「線形重回帰」は、基本的に「まっすぐな関係」、つまり直線や平面しか表現できません。 しかし、実際の現象は「まっすぐな関係」だけでなく、「曲がった関係」を導入しないと説[…]
応答曲面法のPythonでの実装準備
応答曲面モデルは線形重回帰をベースにしており、説明変数の主効果だけでなく、多項式項や交互作用項を加えたモデルです。
説明変数が３つ、多項式項が二次項までの時の応答曲面モデルの式は以下の通りです。
化学で応答曲面法を取り扱う場合、多項式項は多くて三次項までという認識を筆者はしています。
それは、化学において四次以上の複雑な現象はそう多くはないと考えているためです。
まずは、このような応答曲面モデルによるフィッティングをPythonで実装していきましょう。
サンプルデータのダウンロードと想定事例
まずはサンプルデータを以下からダウンロードしてください。
このデータは、とある化学反応を想定して作ったトイデータです。
２つの基質と触媒を溶媒中で混合して目的物を得る化学反応で、副生物も生成します。
そのため、目的物の収率が最も高くなる条件を探索したい、というモチベーションを想定してください。
ここで、説明変数は「基質濃度」「反応温度」「触媒量」、目的変数は「目的物の収率」として解析を進めていきます。
トイデータもこの４つの変数が各列に格納されており、データ数は30となっています。
Pythonによる回帰をやってみよう
それでは、いよいよPythonによる応答曲面法での回帰を行っていきましょう。
Pythonによる基本的な線形回帰は以下記事で詳説していますので、応答曲面法に独特な部分以外は簡単な解説に留めます。
線形回帰はたった2行で実装できる 化学研究者が最も一般的に用いるモデリング手法は何を隠そう「線形回帰」です。 本記事では、この線形回帰をPythonでどのように実装するかを解説します。 実際に手を動かしながら学習できるようにサンプ[…]
（※本記事のサンプルコードではimportを都度行います。これはどの作業でどのライブラリが使われているかを明示するためです。実際はプログラムの最初にまとめて行うのがベター！）
データ読み込み
まずはデータをpandasのデータフレームで読み込みます。
import pandas as pd
folder_path = "データファイルを保管しているフォルダパス" # 適宜変更
# csvファイルをデータフレームで読み込み
df = pd.read_csv(folder_path + "/" + "RSM_サンプルデータ.csv")主効果の中心化
応答曲面法によるモデリングを行うにあたって、説明変数として多項式項と交互作用項（これらをまとめて”高次項”と表現します）を用意します。
その前に、主効果の中心化を行うことが重要です。
多項式項は説明変数の二乗、交互作用は２つの説明変数同士の積、となりますが、主効果をそのままの状態で高次項を算出すると、当然ですが主効果との相関が非常に強くなります。
説明変数同士の相関が強いということは、多重共線性の問題につながります。
そこで、主効果に中心化（平均を差し引いて、平均を0にする）という操作を施してから、高次項を算出します。
そうすることで主効果との相関を消すことができます。
この中心化の効果は別の記事で詳説予定です。
では、中心化の具体的なPythonコードを記述します。
# 目的変数の列名
target_column = "yield(%)"
# 説明変数の列名リストを自動取得（目的変数以外）
explanatory_vars = df.drop(columns=target_column).columns.tolist()
# 中心化用にデータフレームをコピーで新たに作成
df_centering = df.copy()
# 主効果を中心化（各データから平均値を引く）
for var in explanatory_vars:
df_centering[var] = df_centering[var] - df_centering[var].mean()上記コードで「df_centering」という中心化したデータフレームを新たに作成しています。
describe()関数で中心化の確認をしてみましょう。
print(df.describe())
print(df_centering.describe())
df_centeringでは平均がほぼ0になっていることが確認できましたね。
多項式項と交互作用項の導入
では、中心化した主効果を使って高次項を算出していきます。
ここは単純な四則演算で、高次項を新たな列としてデータフレームに追加します。
# 二次項の追加
for var in explanatory_vars:
df_centering[f"{var}^2"] = df_centering[var] ** 2
# 交互作用項の追加
import itertools
for var1, var2 in itertools.combinations(explanatory_vars, 2):
df_centering[f"{var1}×{var2}"] = df_centering[var1] * df_centering[var2]
# 結果確認
print(df_centering.head())交互作用項は説明変数の組み合わせになるので、itertoolsライブラリのcombinations()関数を用いて追加しています。
最後のprint()関数で高次式の列が追加されているのが確認できます。
回帰
ここまでくれば、あとは通常の線形重回帰と流れは同じです。
まず、特徴量行列と目的変数行列に変換して、標準化をしてから回帰を行います。
この流れの詳説は以下の記事を参考にしてください。
線形回帰はたった2行で実装できる 化学研究者が最も一般的に用いるモデリング手法は何を隠そう「線形回帰」です。 本記事では、この線形回帰をPythonでどのように実装するかを解説します。 実際に手を動かしながら学習できるようにサンプ[…]
実際のコードは以下の通りです。
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler
# 説明変数と目的変数に分ける
X = df_centering.drop(columns=target_column).to_numpy()
y = (
df_centering[target_column].to_numpy().reshape(-1, 1)
) # reshape必須（1次元 -> 2次元）
# スケーラーの準備と標準化
X_scaler = StandardScaler()
y_scaler = StandardScaler()
X_scaled = X_scaler.fit_transform(X)
y_scaled = y_scaler.fit_transform(y)
# 回帰
model = LinearRegression()
model.fit(X_scaled, y_scaled)結果の確認
ここも通常の線形重回帰と同じになります。
ここでは精度指標の算出と、「実測値 vs 予測値」の2次元プロットの確認を行っていきましょう。
まず、決定係数とRSMEを下記コードで算出してみましょう。
from sklearn.metrics import mean_squared_error, r2_score
# 予測値の取得
y_pred_scaled = model.predict(X_scaled)
# 元のスケールに戻す
y_pred = y_scaler.inverse_transform(y_pred_scaled)
# R²とRMSEの算出
r2 = r2_score(y, y_pred)
rmse = np.sqrt(mean_squared_error(y, y_pred))
print(f"決定係数 R²: {r2:.4f}")
print(f"RMSE: {rmse:.4f}")決定係数が0.87程度で、比較的良い当てはまりでフィッティングできていることが確認できます。
次に、「実測値 vs 予測値」の2次元プロットを以下コードで描いてみましょう。
import matplotlib.pyplot as plt
# 散布図の描画
plt.figure(figsize=(6, 6)) # 正方形のプロット
plt.scatter(y, y_pred, color="blue", alpha=0.7)
plt.plot([y.min(), y.max()], [y.min(), y.max()], "r--") # 対角線
plt.xlabel("Actual", fontsize=12)
plt.ylabel("Prediction", fontsize=12)
plt.title("Prediction vs Acutual", fontsize=14)
plt.grid(True)
plt.axis("square") # 正方形スケールに調整
plt.tight_layout()
plt.show()以下のようなグラフが得られ、決定係数から示唆されたように比較的対角線上にプロットが乗っており、精度よくデータを予測できていることが分かります。
応答曲面モデルを解釈しよう
ここまでで、交互作用と二次項を含む応答曲面モデルでの線形回帰をPythonで実装しました。
そして、比較的精度の高いモデルを作成することができました。
化学分野でのデータサイエンス活用においては、精度の高いモデルを組んで終わりではなく、そのモデルを正しく解釈して、対象とする現象の知見を得ることが非常に重要です。
そうすることで、さらに良い条件の探索、さらに性能の良い材料開発、といった次のステップにつながります。
応答曲面法では交互作用や多項式といった特殊な変数を含んでおり、通常の線形重回帰よりも解釈が難解なので丁寧に解説していきます。
ということで、作成したモデルを深堀していきましょう。
各説明変数の目的変数への影響度
モデルを解釈するために興味があるのは、「各説明変数がどの程度目的変数に影響するか」「どの説明変数を変えれば目的変数を大きく変えられるか」といったところですよね。
これは、標準化した状態での各説明変数の係数の大小を比較することで確認できます。
標準化状態、つまり説明変数のスケールが統一された状態での係数を比較することで、目的変数への影響が大きい説明変数が何かを考察できます。
モデルの係数はcoef_属性で簡単にアクセスできます。
今回は係数の大小を比較しやすいように可視化してみましょう。
# モデルの係数（標準化した状態）
coefficients = model.coef_.flatten() # flatten()関数で1次元化
# 可視化
# カラム名（説明変数の列）を取得
feature_names = df_centering.drop(columns=target_column).columns.tolist()
plt.figure(figsize=(10, 6))
bars = plt.bar(feature_names, coefficients_flat, color="steelblue")
plt.axhline(0, color="gray", linestyle="--")
plt.ylabel("Standardized Coefficient")
plt.title("Relative Importance of Explanatory Variables")
plt.xticks(rotation=45, ha="right")
plt.tight_layout()
# 数値を棒グラフ上に表示
for bar in bars:
height = bar.get_height()
plt.text(
bar.get_x() + bar.get_width() / 2.0,
height,
f"{height:.2f}",
ha="center",
va="bottom",
)
plt.show()
グラフから、主効果に加えて基質濃度の2次項、温度と触媒量の交互作用項の係数が比較的大きいことが示唆されました。
各説明変数の確からしさ
「どの説明変数が目的変数に影響していそうか」を考えるとき、係数の大小比較に加えて、各説明変数の「確からしさ」という視点も重要です。
例えば今回のケースでは、上のグラフから基質濃度と触媒濃度の交互作用項も係数が0.10となっており、大きくはないですが目的変数に影響する可能性を示唆しています。
この交互作用は本当に目的変数に効くのでしょうか？
言い換えると、この交互作用が目的変数に効くという仮説の「確からしさ」はどの程度なのでしょうか？
この「確からしさ」を定量的に指標化したものがp値というものです。
p値の背景には仮説検定という統計学の確率論の考え方が存在しますが、本記事では深く立ち入りません。
p値は0から1までの範囲を取りうる値で、端的に表現すると「ある説明変数がもし目的変数に影響しないと仮定したとき、算出された係数を偶然取りうる確率」です。
つまりp値が非常に小さい時、その係数値は偶然ではなく、本当に目的変数に影響していると考えることができます。
「p値が非常に小さい時」というのを0.05以下と設定するケースが多いです。
ここまで、つらつらと「確からしさ」の概念を説明してきましたが、つまるところ、「どの説明変数が目的変数に影響していそうか」を考えるとき、係数値が比較的大きく、かつ、p値が小さい（例えば0.05以下）変数を探せばOKです。
では、今回のケースでのp値を確認していきましょう。
以下がコードになります。p値も可視化しましょう。
p値の算出のために、statsmodelsライブラリを使って再度回帰を行っています。
import statsmodels.api as sm
# モデルの切片を追加
X_final = sm.add_constant(X_scaled)
# OLS(最小二乗法)モデルの構築とフィッティング
model_sm = sm.OLS(y_scaled, X_final)
results = model_sm.fit()
# p値の取得
p_values = results.pvalues
# 可視化
feature_names.insert(0, "intercept")
plt.figure(figsize=(10, 6))
bars = plt.bar(feature_names, p_values, color="steelblue")
plt.axhline(0, color="gray", linestyle="--")
plt.ylabel("p values")
plt.title("p values of Explanatory Variables")
plt.xticks(rotation=45, ha="right")
plt.tight_layout()
# 数値を棒グラフ上に表示
for bar in bars:
height = bar.get_height()
plt.text(
bar.get_x() + bar.get_width() / 2.0,
height,
f"{height:.2f}",
ha="center",
va="bottom",
)
plt.show()
p値が0.05以下を優位水準とすると、先ほどの基質濃度と触媒濃度の交互作用項は0.3なので、 “偶然” に係数値が比較的高い値になっただけと考えられます。
ちなみに、一番左端のinterceptは切片のことです。
標準化したデータの回帰では、数学的に切片が0になることが保証されているため、p値も必ず1になります。
さて、係数の大小、およびp値から、主効果と基質濃度の2次項、温度と触媒量の交互作用項の影響が大きいことが示唆されました。
この2次項と交互作用項は何を意味するのか、一つずつ詳しく見ていきます。
多項式項
こちらは比較的分かりやすいです。
2次の場合、2次関数的にその説明変数が効いてくることになるので、係数が正の場合は下に凸、負の場合は上に凸の効き方になります。
今回、基質濃度の2次項の係数が負に比較的大きいので、収率に対して上に凸の効き方になっています。
つまり、収率に対して極大値を持つ変数になっていると考えられます。
では、実際に基質濃度の効き方を2次元プロットで確認してみましょう。
コードは以下の通りです。
# 散布図の描画
plt.figure(figsize=(6, 6)) # 正方形のプロット
plt.scatter(df["Conc.(%)"], df["yield(%)"], color="blue", alpha=0.7)
plt.xlabel("Conc.(%)", fontsize=12)
plt.ylabel("Yield(%)", fontsize=12)
plt.ylim(60, 90)
plt.tight_layout()
plt.show()
グラフを見ると、極大を持っているようにも、40%以上の領域で頭打ちになっているようにも見えますね。
このように頭打ちに見える場合でも、応答曲面法ではモデルは2次項として表現しようとします。
頭打ちの傾向を表現する項（対数項等）が含まれないので当然ですね。
そのため、多項式項の影響が大きいと考えられる場合、その変数と目的変数の2次元プロットを描画して、実験している領域で本当に極大・極小値があるのかを確認することが重要です。
交互作用項
こちらは結構難解です。
説明変数同士の積、というと簡単に聞こえますが、少し見方を変えると「片方の説明変数の係数がもう片方の説明変数の取る値によって変わる」とも言えます。
数式で表現すると以下の通りです。
つまり、片方の変数の目的変数への効き方が、もう片方の変数の値によって変わってしまうのが交互作用です。
どのように目的変数への効き方が変わるのか、というのには主に以下４つのパターンがあります。
片方の変数が目的変数にプラスに効いているとき、もう片方の変数の効き方が「より強くなる=相乗効果」「弱くなる=拮抗効果」「無くなる=飽和効果」「逆転する=逆転効果」というパターンに分かれます。
交互作用項の存在が示唆されたとき、交互作用プロットというものを描画して、上のどのパターンの効果なのかを確認することがモデル解釈には重要です。
交互作用プロットは、x軸に片方の変数、y軸に目的変数をとり、もう片方の変数の各レベルを異なる線でプロットしたものです。
文章では分かりづらいので、実際に描画してみましょう。
今回、温度と触媒量の交互作用項の存在が示唆されたのでこちらの交互作用プロットを描画しましょう。
Pythonで交互作用プロットを描画する関数を以下の通り作りました。
長いのでコピペして使ってください。
import math
def plot_interaction(
df: pd.DataFrame,
categorical_col: str,
continuous_col: str,
target_col: str,
n_levels: int = 4,
):
# 指定された変数を、データ数で均等にn_levels個のグループに分ける
# retbins=Trueで、各グループの境界値を取得する
df_copy = df.copy()
df_copy["category"], bins = pd.qcut(
df_copy[categorical_col],
q=n_levels,
labels=False,
retbins=True,
)
# 全プロットで共通の軸範囲を計算
# 描画範囲に少しマージンを持たせる
x_margin = (df[continuous_col].max() - df[continuous_col].min()) * 0.05
y_margin = (df[target_col].max() - df[target_col].min()) * 0.05
global_xlim = (
df[continuous_col].min() - x_margin,
df[continuous_col].max() + x_margin,
)
global_ylim = (df[target_col].min() - y_margin, df[target_col].max() + y_margin)
# プロットのグリッドを計算
# 1行に最大4つまでプロットを配置
n_cols = min(n_levels, 4)
n_rows = math.ceil(n_levels / n_cols)
fig, axes = plt.subplots(n_rows, n_cols, figsize=(n_cols * 4, n_rows * 3.5))
# axesを1次元配列に変換して扱いやすくする
# n_levelsが1の場合にaxesがndarrayにならないため、条件分岐で対応
if n_levels == 1:
axes = np.array([axes])
else:
axes = axes.flatten()
# 各水準（カテゴリ）ごとにプロットを作成
for i in range(n_levels):
ax = axes[i]
# 現在のカテゴリのデータを抽出
subset_df = df_copy[df_copy["category"] == i]
if subset_df.empty:
ax.text(0.5, 0.5, "No data", ha="center", va="center")
ax.set_title(f"{categorical_col} Level {i+1}")
continue
# 散布図をプロット
x_data = subset_df[continuous_col]
y_data = subset_df[target_col]
ax.scatter(x_data, y_data, alpha=0.6)
# 回帰直線を計算してプロット
# データが2点以上ある場合のみ回帰線を引く
if len(x_data) >= 2:
slope, intercept = np.polyfit(x_data, y_data, 1)
# 回帰線用のX値を作成（プロット範囲全体に線を引くため）
reg_x = np.array(global_xlim)
reg_y = slope * reg_x + intercept
ax.plot(reg_x, reg_y, color="red", linestyle="--", linewidth=2)
# 水準の定義をタイトルとして表示
lower_bound = bins[i]
upper_bound = bins[i + 1]
ax.set_title(f"{categorical_col}: {lower_bound:.2f} - {upper_bound:.2f}")
# 軸ラベル
ax.set_xlabel(continuous_col)
ax.set_ylabel(target_col)
ax.grid(True, linestyle="--", alpha=0.5)
# 全プロットで軸範囲を統一
ax.set_xlim(global_xlim)
ax.set_ylim(global_ylim)
# プロットエリアを正方形にする
x_range = global_xlim[1] - global_xlim[0]
y_range = global_ylim[1] - global_ylim[0]
if x_range > 0 and y_range > 0:
ax.set_aspect(x_range / y_range, adjustable="box")
# 使われなかった余分なプロットを非表示にする
for j in range(i + 1, len(axes)):
fig.delaxes(axes[j])
# 全体のタイトル
fig.suptitle(f"Interaction plot", fontsize=16, y=0.98)
# レイアウトを調整
plt.subplots_adjust(wspace=0.4, hspace=0.6, top=0.8, bottom=0.08, left=0.1, right=0.9)
plt.show()この関数は、データセットのデータフレームと、交互作用プロットを描画したい2つの説明変数のカラム名、目的変数のカラム名、そして片方の変数のレベル数を引数に取ります。
最初に渡した説明変数のカラム名のレベルごとのプロットを描画します。
では、実際にこの関数を使って交互作用プロットを描画してみましょう。
ここではレベル数を4つにしてみます。
plot_interaction(df, "Temp.(℃)", "Cat.(mol%)", "yield(%)", 4)
交互作用プロットを見ると、温度領域によって触媒濃度の効き方が変わっていることが確認できますね。
温度が低い時（65.00-72.25℃）は触媒濃度を上げるほど収率が上がるのが明白ですが、温度を上げるとその効果が頭打ちになってきて、温度が高い時（88.75-99.00）は逆に触媒濃度を上げると収率が下がっている傾向が見えています。
つまり今回のケースでは「飽和効果」あるいは「逆転効果」と言えるでしょう。
結論
では、最後にモデルの数式を求めてみましょう。
先ほどは標準化した状態での係数の大小を比較しましたが、係数を元のスケールに戻してみましょう。
コードは以下の通りです。
# 元のスケールでの係数の表示
# Scalerから元の平均と標準偏差を取得
X_mean = X_scaler.mean_
X_scale = X_scaler.scale_
y_mean = y_scaler.mean_[0]
y_scale = y_scaler.scale_[0]
# 元のスケールでの係数計算
coef_original = coefficients * (y_scale / X_scale)
print(coef_original)
# 元のスケールでの切片計算
intercept = model.intercept_[0]
interc_original = y_mean - np.sum(coef_original * X_mean) + intercept
print(interc_original)スケールを戻すための計算は、標準化が「平均を引いて標準偏差で割る」というのを考慮すれば簡単に導出できます。
興味があれば実際に書いて確かめてみてください。
出力された係数から、上述の通りp値も考慮して「主効果」「基質濃度の二次項」「反応温度と触媒量の交互作用項」をモデルに組み込むとすると、数式は以下の通りになります。
ここで、今回求まった係数値はあくまで各変数を中心化した時のものなので、数式も各変数の平均値が入っています。
この数式を使って、x1, x2, x3のある特定の組み合わせの条件での収率を予測することができます。
さらに、収率が最大になるような条件を逆に求めることも可能になります。
まとめ
応答曲面法による線形回帰について、Pythonで実装するための手順を詳しく見てきました。
説明変数を中心化したのちに高次項（多項式項、交互作用項）を追加して、標準化したのちにフィッティングし、モデルの係数の大小やp値から目的変数への影響度を判定しました。
そして、2次項や交互作用項の存在が示唆された場合は、2次元プロットや交互作用プロットを描画して、実際にどのように目的変数に影響しているかを考察しました。
この記事を読むことで、応用的な線形回帰の流れ全体を掴むことができたと思います。
かつ、具体的なPythonコードを記載しているので、すぐにでもご自身の研究に活用できるかと思います！