線形回帰はたった2行で実装できる
化学研究者が最も一般的に用いるモデリング手法は何を隠そう「線形回帰」です。
本記事では、この線形回帰をPythonでどのように実装するかを解説します。
実際に手を動かしながら学習できるようにサンプルデータも用意しています!
scikit-learnというライブラリを使えば、回帰の部分はたった2行で実装できます。
簡単なのでご自身の環境で試しながら読み進めていただくと理解も深まります!
「そもそも線形回帰って何?」「線形回帰は聞いたことがあるけど、ちゃんとは理解できていない」という方は以下の記事をご覧になってください。
「線形回帰」とは?~線形単回帰と線形重回帰~ 実験計画法を用いたデータ駆動型の研究開発は、企業研究者が最短で目標達成するのに必要なスキルです。 データ駆動型の研究プロセスの中で、花形とも言える「モデリング」には様々な手法があります。 […]
また、必要なライブラリをまだインストールしていない方は以下記事を参考にしてください。
一括インストールから用途別に徹底解説 はじめに 化学研究や実験データ解析、さらには自動化や機械学習による開発加速が求められる中、Pythonはこうした業務を大幅に効率化するための強力なツールです。 本記事では、化学[…]
Pythonによる線形回帰の全体像
Pythonによる処理の全体像を図示します。
まずデータを読み込んで、これを回帰するのに適した形に整形します。
そして、実際の回帰(=フィッティング)を行って、最後にモデルの精度を確認します。
本記事では結果の可視化を含めて、解説していきます。
各段階で適したライブラリを使って処理していくことになります。
コードを交えて順に詳しく見ていきます。
データ準備
化学系の場合、データはcsvやエクセルで保管されていることが多いですよね。
これをまずはPythonで扱えるように読み込む必要があります。
csvやエクセルのデータ読み込みはpandasというライブラリが適しています。
pandasで読み込めば、データフレームという表形式のデータ型で扱えるので、その後の処理が楽になります。
pandasによるデータ読み込みの詳しい内容は以下記事が参考になります。
実際にデータを読み込んでみよう
以下ボタンからサンプルデータをダウンロードしてみてください。
ご自身のデータがある場合はそちらを使ってください。
データファイルを適当な場所に保存していただき、以下コードで読み込んでみましょう。
read_csv()関数を使ってデータフレームで読み込みます。
import pandas as pd # pandasのインポート
folder_path = "データファイルを保存したフォルダのパス" # 適宜変更
# csvファイルをデータフレームで読み込み
df = pd.read_csv(folder_path + "/" + "regression_data.csv")
# データフレームの最初の5行を表示(確認作業)
df.head()
こんな感じで表形式でデータが読み込めたことを確認できるはずです。
(※本記事のサンプルコードではimportを都度行います。これはどの作業でどのライブラリが使われているかを明示するためです。実際はプログラムの最初にまとめて行うのがベター!)
データ前処理
NumPy配列への変換
後の工程を考慮すると、説明変数と目的変数を分けて、かつ、それぞれNumPy配列、つまり行列に変換しておいた方が良いです。
名前の通り、numpyを用いて変換できます。変換イメージは下図の通り。
図のように、複数の説明変数を列方向に結合したものを特徴量行列と呼びます。
なぜ行列を用いるのかについて簡単に触れます。
それは単純に処理が早くなるからです。
データの数が多く、かつ変数も増えていくと計算量は肥大化しますが、大量のデータを扱うときにforループなんかを回すとめちゃくちゃ時間がかかります。
そこで、大量のデータを行列形式で格納すれば、一気に計算でき、処理スピードが全然違います。
標準化
次に変数のスケールを合わせる処理を行います。
これは説明変数同士の大きさが異なる場合に、最適化プロセスが不安定になったり、モデルの性能が低下するのを防ぎます。
例えば、変数x1は0~10の間の値で、変数x2は1万~100万の間の値という場合、変数x2の値にモデルが引きずられてしまうことがあります。
なんとなく想像できますね。
そこで、各変数に「標準化(=スケーリング)」という処理を施して、平均が0で分散が1になるように変換することで、変数を平等に評価できるようにします。
目的変数の標準化は線形回帰や他のモデルでは不要なこともありますが、ガウス過程回帰などでは必須になるため、標準化しておく癖をつけましょう。
実際に前処理してみよう
先ほどのデータフレームを使って前処理してみましょう。
まずはNumPy配列への変換です。データフレームのto_numpy()メソッドを使います。
以下コードを試してみましょう。
import numpy as np # numpyのインポート
# 特徴量行列の生成(変数が2つ以上の場合は2重カッコ使用)
X = df[["x1", "x2"]].to_numpy()
print(type(X)) # 型を表示(確認作業)
# 目的変数行列の生成(変数が1つ場合はreshape(-1, 1)でn行 x 1列行列に変換)
y = df["y"].to_numpy().reshape(-1, 1)
print(type(y)) # 型を表示(確認用)
変数が1つのときと、2つ以上の時でデータフレームを行列化するコードが若干異なることに注意してください。
2つ以上の時はカラム指定の際に二重カッコを用い、1つの場合は最後にreshape(-1, 1)をつけてください。
気になる方はX, yもprint関数で確認してみましょう。
次に標準化を行います。
標準化にはScikit-learnのStandardScalerのfit_transform()メソッドを用います。
# StandardScalerのインポート
from sklearn.preprocessing import StandardScaler
# Xの標準化
scaler_X = StandardScaler() # インスタンス化
X_scaled = scaler_X.fit_transform(X) # Xのスケーリング実行
# yの標準化
scaler_y = StandardScaler()
y_scaled = scaler_y.fit_transform(y) # yのスケーリング実行
# スケーリングの結果確認(確認作業)
print("スケーリング後のX平均:", np.mean(X_scaled, axis=0))
print("スケーリング後のX分散:", np.var(X_scaled, axis=0))
print("スケーリング後のy平均:", np.mean(y_scaled, axis=0))
print("スケーリング後のy分散:", np.var(y_scaled, axis=0))
スケーリング後の各説明変数の平均と分散を確認すると、ほぼ0と1になっているのが確認できます。
回帰
ここまでくれば回帰は一瞬です。
標準化した特徴量行列と目的変数行列をScikit-learnのLinearRegressionオブジェクトのfitメソッドに引数としてを渡すだけです。
fitメソッドの中で最小二乗法によるフィッティングが行われ、各説明変数の係数と切片の値がLinearRegressionオブジェクトの属性として格納されます。
実際に回帰してみよう
回帰はたったの2行です!
これで最小二乗法によるフィッティングが完了します。
# Scikit-learnのLinearRegressionクラスをインポート
from sklearn.linear_model import LinearRegression
model = LinearRegression() # インスタンス化
model.fit(X_scaled, y_scaled) # 回帰の実行
モデルの確認
たった2行でフィッティングしたモデルが実際にどんなものか見ていきましょう。
ここでは主に2つの確認方法について解説します。
本来、作ったモデルは様々な側面で評価しますが、今回は基本的なところに留めます。
モデル精度の確認
ここでいうモデルの精度は「実測値とモデル予測値の一致性」のことです。
一致性には様々な指標がありますが、R²値(決定係数)やRMSE(二乗平均平方根誤差)が代表的です。
また、視覚的には「実測値 vs 予測値」の2次元プロットもよく用いられます。
予測値を縦軸、実測値を横軸に取り、プロットが対角線上に位置しているほど一致性が高いことが読み取れます。
係数、切片の確認
精度の次に、フィッティングにより算出された係数と切片を確認します。
この係数から、各説明変数がどのように目的変数に影響しているかを考察できます。
「この説明変数は目的変数にほとんど影響しない」「この説明変数をこう変えれば目標に到達できそう」といったことを係数から読み取ります。
実際にモデルを確認してみよう
モデルの評価を行うために、まずはpredict()メソッドを使用して、作成したモデルの予測値を出力します。
この時、標準化の逆の操作、つまり元のスケールに戻すのを忘れないようにしてください。
標準化した時のStandardScalerオブジェクトのinverse_transform()メソッドをを使用します。
R²値とRMSEの算出は、Scikit-learnのr2_score()とmean_squared_error()を用います。
from sklearn.metrics import mean_squared_error, r2_score
# 予測値の取得
y_pred_scaled = model.predict(X_scaled)
# 元のスケールに戻す
y_pred = scaler_y.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}")
mean_squared_error()では平均二乗誤差(MSE)を算出する関数なので、sqrt()関数で平方根を取ってRMSEにしています。
次に、「実測値 vs 予測値」の2次元プロットを描いてみましょう。
matplotlibというライブラリを使って、以下コードを試してみてください。
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) # x軸ラベル設定
plt.ylabel("Prediction", fontsize=12) # y軸ラベル設定
plt.title("Prediction vs Acutual", fontsize=14) # グラフタイトル設定
plt.grid(True) # グリッド線の表示
plt.axis("square") # 正方形スケールに調整
plt.tight_layout()
plt.show()
このように2次元プロットが表示されるはずです。
対角線上にプロットが乗っているほど実測データをうまく予測できていることになります。
今回はR²値が0.85程度なので、比較的きれいに対角線上に乗っていますね。
最後に係数と切片を算出してみましょう。
以下コードで標準化した状態での数値と、標準化を元に戻した時の数値を表示できます。
係数はモデルのcoef_属性、切片はintercept属性にそれぞれ格納されています。
# 標準化したままの説明変数の係数
coefficients = model.coef_
print(coefficients)
# 標準化したままの切片
intercept = model.intercept_[0]
print(intercept)
# 元のスケールでの係数の表示
# Scalerから元の平均と標準偏差を取得
X_mean = scaler_X.mean_
X_scale = scaler_X.scale_
y_mean = scaler_y.mean_[0]
y_scale = scaler_y.scale_[0]
# 元のスケールでの係数計算
coef_original = coefficients * (y_scale / X_scale)
print(coef_original)
# 元のスケールでの切片計算
interc_original = y_mean - np.sum(coef_original * X_mean)
print(interc_original)
元のスケールでの係数や切片は、標準化前の平均や標準偏差を使って計算されます。
計算式の説明は割愛します。
元のスケールでの係数はx1, x2ともに0.33程度、切片は1程度になりました。
実はこのデータセットは以下数式を元に、適当に誤差を発生させて作ったものです。
線形回帰によって比較的精度よく予測できていることが確認できますね。
おまけ:3Dプロットによる可視化
おまけで今回作成したモデルが3次元空間でどのように見えるかを可視化するコードを紹介します。
説明変数が2つのため、plotlyというライブラリを用いて、以下コードで3次元空間を表現できます。
マウスで3Dグラフをぐりぐり動かせるので、モデルが作る平面が実験点を比較的うまくフィッティングできていることを確認してみてください。
import plotly.graph_objects as go
# 回帰面の作成範囲(x1, x2 は元のスケールで指定)
x1_grid = np.linspace(df["x1"].min(), df["x1"].max(), 30)
x2_grid = np.linspace(df["x2"].min(), df["x2"].max(), 30)
x1_mesh, x2_mesh = np.meshgrid(x1_grid, x2_grid)
# 予測用グリッド点を作成し、スケーリング
X_pred_original = np.column_stack((x1_mesh.ravel(), x2_mesh.ravel()))
X_pred_scaled = scaler_X.transform(X_pred_original)
# 予測(スケーリングされた y_pred を元スケールに戻す)
y_pred_scaled = model.predict(X_pred_scaled)
y_pred_original = scaler_y.inverse_transform(y_pred_scaled).reshape(x1_mesh.shape)
# プロット
fig = go.Figure()
# データ点(元のスケール)の3D散布図
fig.add_trace(
go.Scatter3d(
x=df["x1"],
y=df["x2"],
z=df["y"],
mode="markers",
marker=dict(size=2, color="#1e2e53"),
name="Data",
)
)
# 回帰面の描画(元スケール)
fig.add_trace(
go.Surface(
x=x1_mesh,
y=x2_mesh,
z=y_pred_original,
opacity=0.7,
showscale=False,
colorscale=[[0, "blue"], [1, "blue"]],
name="Regression Surface",
)
)
# 軸範囲はデータに合わせて自動設定も可能ですが、固定したい場合は以下を調整
fig.update_layout(
title="Multiple Linear Regression: Interactive 3D Surface (Original Scale)",
scene=dict(
xaxis_title="x1",
yaxis_title="x2",
zaxis_title="y",
xaxis=dict(range=[df["x1"].min(), df["x1"].max()]),
yaxis=dict(range=[df["x2"].min(), df["x2"].max()]),
zaxis=dict(range=[df["y"].min(), df["y"].max()]),
),
margin=dict(l=0, r=0, t=40, b=0),
width=900,
height=700,
)
fig.show()
おわりに
今回は線形回帰をPythonで実装する方法を具体的なコードでお示ししました。
データを読み込んで、適切に変形(行列変換、標準化)をし、回帰を行って、最後に結果を確認するところまで丁寧に解説しました。
結果の確認ではグラフを使った可視化が非常に重要です。
可視化することで初めて見えてくることもあるので、絶対に実施するようにしましょう!
全コード
import pandas as pd
import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score
import matplotlib.pyplot as plt
folder_path = "データファイルを保存したフォルダのパス" # 適宜変更
# csvファイルをデータフレームで読み込み
df = pd.read_csv(folder_path + "/" + "regression_data.csv")
# データフレームの最初の5行を表示(確認作業)
df.head()
# 特徴量行列の生成(変数が2つ以上の場合は2重カッコ使用)
X = df[["x1", "x2"]].to_numpy()
# 目的変数行列の生成(変数が1つ場合はreshape(-1, 1)で列ベクトルに変換)
y = df["y"].to_numpy().reshape(-1, 1)
# Xの標準化
scaler_X = StandardScaler()
X_scaled = scaler_X.fit_transform(X)
# yの標準化
scaler_y = StandardScaler()
y_scaled = scaler_y.fit_transform(y)
print("スケーリング後のX平均:", np.mean(X_scaled, axis=0))
print("スケーリング後のX分散:", np.var(X_scaled, axis=0))
print("スケーリング後のy平均:", np.mean(y_scaled, axis=0))
print("スケーリング後のy分散:", np.var(y_scaled, axis=0))
# モデリング
model = LinearRegression()
model.fit(X_scaled, y_scaled)
# 予測値の取得
y_pred_scaled = model.predict(X_scaled)
# 元のスケールに戻す
y_pred = scaler_y.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}")
# 散布図の描画
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()
# 標準化したままの説明変数の係数
coefficients = model.coef_
print(coefficients)
# 標準化したままの切片
intercept = model.intercept_[0]
print(intercept)
# 元のスケールでの係数の表示
# Scalerから元の平均と標準偏差を取得
X_mean = scaler_X.mean_
X_scale = scaler_X.scale_
y_mean = scaler_y.mean_[0]
y_scale = scaler_y.scale_[0]
# 元のスケールでの係数計算
coef_original = coefficients * (y_scale / X_scale)
print(coef_original)
# 元のスケールでの切片計算
interc_original = y_mean - np.sum(coef_original * X_mean)
print(interc_original)