SPSS ModelerでPLSを実行したい
みなさん、こんにちわ。今回は、部分最小二乗(partial least squares, PLS)回帰をSPSS Modelerで挑戦してみたいと思います。このPLSに関してもユーザーの方からの問い合わせがありました。
Rでの実装方法はすでに紹介されていますので、私は、Pythonで実装してみようと思います。
このPLSの実装については、複数回に分けて紹介する予定です。
今回は、PLSでの潜在変数(PLS成分)の最適な数を調査します。
次回以降に、ModelerでPLSをモデル化しようと思います。
※. 3/4 拡張ノードでのモデル精度確認部分を大幅に修正しました。申し訳ございません。
※ 2025/03/14 記事で使用したストリームをGitHubにアップロードしました。
データについては、以下のリンクより"gasoline.csv"をダウンロードしてください。
1. PLS(partial least squares, PLS)とは?
多重共線性(説明変数同士が強く相関する問題)に対処しながら、目的変数を予測するための回帰手法です。主に、説明変数が多く、互いに強く相関している場合や、サンプル数が変数の数より少ない場合に有効です。
PLSの特徴
①. 次元削減と回帰を同時に行う
PCA(主成分分析)と似ており、説明変数の情報を少数の潜在変数(PLS成分)に圧縮し、それを用いて回帰を行う。
これにより、多重共線性の影響を低減しつつ、予測精度を向上させる。
②. 多重共線性に強い
例えば、化学データ解析(スペクトルデータ)などで、非常に多くの相関する変数がある場合に有効。
③. データ数が少なくても使える
説明変数の数がサンプル数より多い(p >> n 問題)場合でも適用可能。
2. さっそく準備
①. 環境情報
今回実装した環境の情報は以下になります。
ソフトウェア等 | バージョン |
---|---|
OS | Windows 11 Pro |
SPSS Modeler | 18.6 (windows) |
Python | 3.10.7 (デフォルト) |
scikit-learn | 1.6.1 |
numpy | 1.26.4 |
pandas | 2.2.3 |
ここで、注意が必要なのはnumpyです。
今回、scikit-learnを導入する時に、自動でnumpy 2.2.3がインストールされました。
ですが、2.2.3の場合、以下のエラーがでてうまく動作しませんでした。
ValueError: numpy.dtype size changed, may indicate binary incompatibility.
そのため、v1.26.4を再導入しました。
※.このNumpyのバージョン問題は、テキストマイニングの時にも発生しました。
①. Scikit-Learnを導入
PLSを使うために、Scikit-Learnを導入しましょう。
下記コマンドは、SPSS Modeler v18.6でのデフォルトパスへのインストール時の例となります。
C:\Program Files\IBM\SPSS\Modeler\18.6\python_venv\Scripts>python.exe -m pip install Scikit-Learn
導入後には、Numpy v1.26.4 を再導入しましょう。
②. データ
以下のデータを使用させていただきました。
NIR スペクトル(近赤外分光データ)から、オクタン価(Octane Number)を予測するのに利用できるデータですね。
オクタン価(Octane Number) は、ガソリンの品質を示す指標の1つ。
数値が高いほど「ノッキング(異常燃焼)」が起きにくい → 高性能なエンジンに適している。
NIR スペクトル(近赤外分光データ) は、近赤外線(800~2500nmの波長範囲)を使い、ガソリンの化学組成を分析する方法で「波長ごとの吸光度(log(1/R))」として記録される。
今回は、gasoline.csvを利用しています。
3. グラフでデータ確認
スペクトルとオクタン値の関係を視覚化してみます。
以下のように、横軸にスペクトル、縦軸に吸光度(log(1/R)を配置したグラフを書いてみます。
そしてオクタン値毎に色分けして視覚化しました。Git Hubのサイトでも3D化して視覚化していますね。
①.データを入力
gasoline.csvを可変長ファイルノードで入力します。ファイルタブでは特に変更する部分はありません。
データタブで、ストレージを変更します。
下記の3つのフィールドが「文字列」となっているのを、「実数」に変更します。
データ型タブで、「値の読み込み」をし、インスタンス化します。
②.スペクトルフィールドを縦持ちに
さて、冒頭のグラフイメージの通り、X軸にスペクトル値を配置したいので、フィールド名を値に変換するために、行列入れ替えでスペクトルフィールド名を縦持ちにします。
行列入替方法は、「フィールドからレコードへ」
索引に、「octane」を設定
値に、スペクトルフィールドを全て選択します。
以下のようになります。フィールド名が列に展開されていますよね。
valueは、もともと、スペクトルフィールドに格納されていた、吸光度(log(1/R)です。
③.フィールド名を変更
フィルターノードで先ほど展開したスペクトルフィールド名をwavelengthに変更します。
④.スペクトルの値を変換
"NIR.900 nm" のような値を "900" のように、文字列を削除して数値だけにします。
置換ノードで先ほどフィールド名を変更した「wavelength」フィールドの値を置換します。
以下のようにします。
・フィールド値の先頭から4文字分を削除 = allbutfirst(4, field )
・フィールド値の後ろから3文字分を削除 = allbutlast(3, field )
・最後、整数型に変換 = to_integer( field )
これらを組み合わせて、
to_integer( allbutlast( 3, allbutfirst( 4, wavelength) ) )
としています。
⑤.オクタン値を名義型に変更
オクタン値で色分けしたいので、数値型ではなく、カテゴリ型である「名義型」に変更します。
ここでは、データ型ノードを使用しています。
⑥.線グラフを書く - 散布図ノードで!!
準備が終わったのでグラフを描画しようと思うのですが、イメージしているようなグラフを線グラフノードや時系列グラフノードではうまく描画できないんですよね。
そこで、登場するのが散布図ノードです。
X軸とY軸に数値フィールドを指定できて、カテゴリフィールドを使って色分けすることができます。
さらに、"散布図"ノードなのに、線グラフも描画できるんです。
まずは、Xフィールドに「wavelength」、Yフィールドに「value」そして、オーバーレイの色に「octane」を指定します。
つづいて、オプションタブで「スタイル」を線に変更します。これで線グラフが描画できるんです。
ちなみに、X軸、Y軸の値の範囲を指定できるのも便利なところです。
ぜひ活用してみてください。
グラフを観察すると、オクタン値毎に吸光度(log(1/R)に微妙な差異がありそうですね。
差が出るスペクトルと出ないスペクトルもありそうです。
いよいよ、PLSでオクタン値を予測してみます。
4. PLSモデルを複数作成して最適な潜在変数の数を見つける
さて、いよいよPLSモデルを作成していくのですが、PLSは、説明変数の情報を少数の潜在変数(PLS成分)に圧縮し、それを用いて回帰を行うものです。
この潜在変数の数をいくつにすれば、精度がよくなるのか?が分からないのです。
そこで今回は潜在変数のパラメータ値を変えながらモデルを作成して、スコアの高い潜在変数の値を見つけようと思います。
この工程の順番は
- データ入力
- 401個の説明変数の絞り込み
- PLSモデルの作成及びパラメータ値の調査
とりなます。
①.データを入力
gasoline.csvを可変長ファイルノードで入力します。
ここは、先ほどグラフを描画した時と同じ設定をしてください。
NIR 1246 nm
NIR 1248 nm
NIR 1348 nm
上記3つのフィールドのストレージを"実数"に変更するのをお忘れなく。
②.学習データとテストデータに分割
データ区分ノードで学習データとテストデータに分割します。
ここでは、学習データ 80 : テストデータ 20 と設定します。
③. 不要なフィールドとレコードの削除
条件抽出ノードでデータは学習データのみにします。
フィルターノードでoctaneとスペクトルフィールド以外は削除します。
(以下のとおり、2つだけです。)
④. ロールの指定
データ型ノードでoctaneを「対象」、スペクトルフィールド全てを「入力」とします。
⑤. 説明変数の絞り込み
さて今回のデータは、
レコード数 : 60
スペクトルフィールド数 : 401
というもので、従来の重回帰分析を使おうとしても自由度が全然足りなくて解析できません。
そこで、今回はモデル作成タブにある特徴量選択ノードを使い、octaneと相関の高い変数にまずは絞り込みます。
特に設定は変更せずに実行してください。
するとモデルナゲットができます。内容を確認するとoctane(目的変数)に対して相関の高い変数にチェックが入った状態になります。
閾値に達していない変数は自動的にスクリーニングされます。
※.分析者が自由にチェックを入れたり外したりできます。
今回は61フィールドが選択された状態になりました。
401→61と340ほど説明変数を削減できました。
くわしい使い方は、西牧さんの記事も参考にしてください。
さてここから拡張ノードでPLSを実装していきます。
⑥. 拡張の出力ノードで最適な潜在変数の数を探す
この記事では、モデルを実装するというより、まずは最適なパラメータを探すことが目的なので、
拡張モデルノードではなく、拡張の出力ノードを使います。
ここでは、
- 潜在変数の数を1・2・3と増やしながらモデルを構築
- それぞれのモデルを交差検証して評価
- 一番評価の良かったモデルの潜在変数の数を確認
という流れでプログラムを書いています。
また、説明変数と、目的変数を変更すれば他のデータでも流用可能です。
i. Pythonコード全体
まずは、Pythonのコード全体です。
#-------------------------------------------------------
# ライブラリインポート
#-------------------------------------------------------
# Modeler用ライブラリ
import modelerpy
# NumpyとPandas
import numpy as np
import pandas as pd
# sklearn - PLS
from sklearn.cross_decomposition import PLSRegression
# cross validation
from sklearn.model_selection import cross_val_predict, KFold
# 精度検証用 R2乗、MSE、MAE
from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error
#-------------------------------------------------------
# データ入力
#-------------------------------------------------------
# 事前にデータ区分で学習・テストデータは分けているのでそのまま入力
# Pandasで入力
modelerData = modelerpy.readPandasDataframe()
# 説明変数 - NIR.900 rm以降の列を全て
xdata = modelerData.iloc[:, 1:] # 2列目(index=1)以降を選択
# 目的変数 - octane
ydata = modelerData['octane']
#-------------------------------------------------------
# データ加工
#-------------------------------------------------------
# PLSで使用するために、説明変数をNumpyのndarrayに変換
input_data = xdata.values
# 目的変数はPandas.SeriesなのでそのままでOK
target_data = ydata
#-------------------------------------------------------
# 交差検証の分割数を定義
# ここも値など適宜変更していください。
#-------------------------------------------------------
cnt_cv=5
# KFold で交差検証
kf = KFold(n_splits=cnt_cv, shuffle=True, random_state=42)
#-------------------------------------------------------
# n_componentsの最大数を定義
#-------------------------------------------------------
#レコード数 or 説明変数 の小さい方
max_comp = min(input_data.shape[0], input_data.shape[1])
#交差検証をするため、n_componentsの数の制限が変わる
n_samples_train = int(input_data.shape[0] * ( (cnt_cv - 1) / cnt_cv ) )
max_comp = min( n_samples_train, max_comp )
#精度格納用配列
r2_scores = [] #R2乗用
mse_scores = [] #MSE用
mae_scores = [] #MAE用
#-------------------------------------------------------
# PLS構築開始
# 1. 最初のループ
# 主成分数を1つずつ増やしながらPLSを構築
# 2. 1で構築したモデルでScoreを計算して配列に格納しておく
#
#-------------------------------------------------------
# n_componentsの最大数分(max_comp)ループさせてPLSを構築
# PLSで使用する主成分数を増やしながら複数のモデルを構築
for cnt_ncomp in range(max_comp):
# ループ数をn_componentsに設定
pls_base = PLSRegression( n_components = cnt_ncomp + 1 )
# PLSモデル実行
pls_base.fit( input_data, target_data )
# 交差検証 :分割数は5
#target_pred = cross_val_predict(pls_base, input_data, target_data, cv=cnt_cv)
# 交差検証 :KFold適用
target_pred = cross_val_predict(pls_base, input_data, target_data, cv=kf)
# 精度計算
#R2乗
r2 = r2_score(target_data, target_pred)
#MSE
mse = mean_squared_error(target_data, target_pred)
#MAE
mae = mean_absolute_error(target_data, target_pred)
#配列に格納
r2_scores.append(r2) #R2乗
mse_scores.append(mse) #MSE
mae_scores.append(mae) #MAE
#-------------------------------------------------------
# 精度の確認
# R2乗<最大値>、MSE & MAE <最小値> を特定
#-------------------------------------------------------
# 最適な n_components を決定
# R2乗 最大値
best_n_r2 = np.argmax(r2_scores) + 1 # R^2 が最大の主成分数
best_r2 = max(r2_scores)
# MSE 最小値
best_n_mse = np.argmin(mse_scores) + 1 # MSE が最小の主成分数
best_mse = min(mse_scores)
# MAE 最小値
best_n_mae = np.argmin(mae_scores) + 1 # MAE が最大の主成分数
best_mae = min(mae_scores)
# 結果出力
print(f"R^2 MAX n_components: {best_n_r2}, R^2: {best_r2:.4f}")
print(f"MSE MIN n_components: {best_n_mse}, MSE: {best_mse:.4f}")
print(f"MAE MIN n_components: {best_n_mae}, MAE: {best_mae:.4f}")
ii. Pythonコード詳細
では、詳細をみていきましょう。
a. ライブラリのインポート
まずは、PLSと交差検証用のもの含め、必要なライブラリをインポートしています。
※3/4 大幅に修正しました。申し訳ございません。
交差検証用に KFold
精度確認用に r2_score, mean_squared_error, mean_absolute_error を追加でインポートしています。
#-------------------------------------------------------
# ライブラリインポート
#-------------------------------------------------------
# Modeler用ライブラリ
import modelerpy
# NumpyとPandas
import numpy as np
import pandas as pd
# sklearn - PLS
from sklearn.cross_decomposition import PLSRegression
# cross validation
from sklearn.model_selection import cross_val_predict, KFold
# 精度検証用 R2乗、MSE、MAE
from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error
b. データ入力
データをModelerから入力して、説明変数と目的変数に分けています。
コメントアウトしているコードもありますが、どちらでも可能です。
やりやすい方法で指定しましょう。
#-------------------------------------------------------
# データ入力
#-------------------------------------------------------
# 事前にデータ区分で学習・テストデータは分けているのでそのまま入力
# Pandasで入力
modelerData = modelerpy.readPandasDataframe()
# 説明変数 - NIR.900 rm以降の列を全て
#xdata = modelerData.loc[:, 'NIR.900 nm':]
xdata = modelerData.iloc[:, 1:] # 2列目(index=1)以降を選択
# 目的変数 - octane
ydata = modelerData['octane']
#ydata = modelerData.iloc[:, 0] # 1列目(index=0)を選択
c. モデル投入用データの作成
PLSRegressionへ渡すために、説明変数のデータは値のみにしています。
#-------------------------------------------------------
# データ加工
#-------------------------------------------------------
# PLSで使用するために、説明変数をNumpyのndarrayに変換
input_data = xdata.values
# 目的変数はPandas.SeriesなのでそのままでOK
target_data = ydata
d. 各種パラメータの設定
ここは、ちょっとだけ注意が必要な部分です。
※ 3/4 大幅に修正しました。
・交差検証にKFoldど採用しました。
・精度確認用配列をR2乗、MSE、MAE用にしました。
#-------------------------------------------------------
# 交差検証の分割数を定義
# ここも値など適宜変更していください。
#-------------------------------------------------------
cnt_cv=5
# KFold で交差検証
kf = KFold(n_splits=cnt_cv, shuffle=True, random_state=42)
#-------------------------------------------------------
# n_componentsの最大数を定義
#-------------------------------------------------------
#レコード数 or 説明変数 の小さい方
max_comp = min(input_data.shape[0], input_data.shape[1])
#交差検証をするため、n_componentsの数の制限が変わる
n_samples_train = int(input_data.shape[0] * ( (cnt_cv - 1) / cnt_cv ) )
max_comp = min( n_samples_train, max_comp )
#精度格納用配列
r2_scores = [] #R2乗用
mse_scores = [] #MSE用
mae_scores = [] #MAE用
最初の部分は交差検証のための分割数は 5 に設定しています。
KFold(shuffle=True) はデータの偏りを減らし、より一般化されたモデルを作るのに役立つといわれています。
※. 3/4 KFoldを採用しました。大変申し訳ございません。
#-------------------------------------------------------
# 交差検証の分割数を定義
# ここも値など適宜変更していください。
#-------------------------------------------------------
cnt_cv=5
# KFold で交差検証
kf = KFold(n_splits=cnt_cv, shuffle=True, random_state=42)
次が分かりにくいのですが、
PLSRegressionでは、以下のように潜在変数の数を指定する必要があります。
"n_components = 数値 "
但し、指定できる最大値に制限があり、
- モデルに投入するレコード数以下
- モデルに投入する説明変数以下
となります。
そのため、まず最大値はレコード数 or 説明変数の数のうち小さい方を設定しています。
#-------------------------------------------------------
# n_componentsの最大数を定義
#-------------------------------------------------------
#レコード数 or 説明変数 の小さい方
max_comp = min(input_data.shape[0], input_data.shape[1])
ただ、それだけでは今回の場合エラーとなってしまいます。
なぜなら、交差検証をするからです。
交差検証をする場合は、1 / 分割数 のデータを検証用に取っておくため
モデルに投入されるレコードは 分割数 - 1 / 分割数 となります。
例) 分割数が 5 の場合
・検証用データ = 1 / 5
・モデル投入データ = 4 / 5
となります。
そのため、ここではモデル投入データの数を計算し、レコード数 or 説明変数の数の小さい方とさらに比較してい最大値を設定しているのです。
#交差検証をするため、n_componentsの数の制限が変わる
n_samples_train = int(input_data.shape[0] * ( (cnt_cv - 1) / cnt_cv ) )
max_comp = min( n_samples_train, max_comp )
今回の例ですと
max_comp = min(input_data.shape[0], input_data.shape[1])
は
45(レコード数)と61(説明変数の数)なので、
max_comp = 45
その次に交差検証用のレコード数を計算すると
int(input_data.shape[0] * ( (cnt_cv - 1) / cnt_cv ) ) = int(45 * ( (5 - 1) / 5 ) )
= 36 となるので、
max_comp = 36が設定されます。
そして最後に精度格納用の配列を定義
#精度格納用配列
r2_scores = [] #R2乗用
mse_scores = [] #MSE用
mae_scores = [] #MAE用
e. PLSモデル作成+交差検証
さて、PLSモデルを作成してから、交差検証をし、スコアを保存していきます。
※3/4 大幅に修正しました。申し訳ございません。
・交差検証をKFoldをパラメータ指定するものに修正
・精度確認にをR2乗、MSE、MAEを採用
#-------------------------------------------------------
# PLS構築開始
# 1. 最初のループ
# 主成分数を1つずつ増やしながらPLSを構築
# 2. 1で構築したモデルでScoreを計算して配列に格納しておく
#
#-------------------------------------------------------
# n_componentsの最大数分(max_comp)ループさせてPLSを構築
# PLSで使用する主成分数を増やしながら複数のモデルを構築
for cnt_ncomp in range(max_comp):
# ループ数をn_componentsに設定
pls_base = PLSRegression( n_components = cnt_ncomp + 1 )
# PLSモデル実行
pls_base.fit( input_data, target_data )
# 交差検証 :分割数は5
#target_pred = cross_val_predict(pls_base, input_data, target_data, cv=cnt_cv)
# 交差検証 :KFold適用
target_pred = cross_val_predict(pls_base, input_data, target_data, cv=kf)
# 精度計算
#R2乗
r2 = r2_score(target_data, target_pred)
#MSE
mse = mean_squared_error(target_data, target_pred)
#MAE
mae = mean_absolute_error(target_data, target_pred)
#配列に格納
r2_scores.append(r2) #R2乗
mse_scores.append(mse) #MSE
mae_scores.append(mae) #MAE
まずはループ部分ですが、0から先ほど設定したmax_comp(潜在変数の最大値)までcnt_ncompの値をカウントアップさせています。
# n_componentsの最大数分(max_comp)ループさせてPLSを構築
# PLSで使用する主成分数を増やしながら複数のモデルを構築
for cnt_ncomp in range(max_comp):
次に、n_componentsパラメータに cnt_ncomp + 1(cnt_ncompは0からはじまるため) を指定して、PLSモデルを作成しています。
# ループ数をn_componentsに設定
pls_base = PLSRegression( n_components = cnt_ncomp + 1 )
# PLSモデル実行
pls_base.fit( input_data, target_data )
最後に、cross_val_predictを使ってスコアを算出しています。
※. 3/4 大幅に修正しました。
・KFoldを指定してしました。
・R2乗、MSE、MAEを計算して配列に格納
# 交差検証 :分割数は5
#target_pred = cross_val_predict(pls_base, input_data, target_data, cv=cnt_cv)
# 交差検証 :KFold適用
target_pred = cross_val_predict(pls_base, input_data, target_data, cv=kf)
# 精度計算
#R2乗
r2 = r2_score(target_data, target_pred)
#MSE
mse = mean_squared_error(target_data, target_pred)
#MA
mae = mean_absolute_error(target_data, target_pred)
#配列に格納
r2_scores.append(r2) #R2乗
mse_scores.append(mse) #MSE
mae_scores.append(mae) #MAE
f. 精度の一番良い n_components を出力
※. 3/4 大幅に修正しました。申し訳ございません。
最後は、各種精度が格納された配列から一番スコア高いものを探して終わりです。
np.argmax(配列) は、配列の値で一番大きい値のインデックスを返してくれます。
np.argmin(配列) は、配列の値で一番小さい値のインデックスを返してくれます。
R2乗 Scores → 各分割ごとの 決定係数のスコア(高い方がよい)
MSE / MAE → 誤差指標の平均値(小さいほど良い)
#-------------------------------------------------------
# 精度の確認
# R2乗<最大値>、MSE & MAE <最小値> を特定
#-------------------------------------------------------
# 最適な n_components を決定
# R2乗 最大値
best_n_r2 = np.argmax(r2_scores) + 1 # R^2 が最大の主成分数
best_r2 = max(r2_scores)
# MSE 最小値
best_n_mse = np.argmin(mse_scores) + 1 # MSE が最小の主成分数
best_mse = min(mse_scores)
# MAE 最小値
best_n_mae = np.argmin(mae_scores) + 1 # MAE が最大の主成分数
best_mae = min(mae_scores)
# 結果出力
print(f"R^2 MAX n_components: {best_n_r2}, R^2: {best_r2:.4f}")
print(f"MSE MIN n_components: {best_n_mse}, MSE: {best_mse:.4f}")
print(f"MAE MIN n_components: {best_n_mae}, MAE: {best_mae:.4f}")
iii. 結果の確認
※ 3/4 大幅に修正しました。申し訳ございません。
拡張の出力ノードを実行すると、以下のように出力されます。
n_componentsの最適な数は、6 のようですね。
最終的に、説明変数の数は
- データ入力時 401
- 特徴量選択字 61
- PLSで潜在変数化 6
といったように、最終的に401 → 6 となりました。
ちなみに、KFoldを使わない場合で、分割数が 5 の場合は、以下のように
n_componentsは 7 で若干ですが、各指標の値が低くなりました。
5. まとめ
今回は、ユーザーの方より質問のあった、PLSモデルに挑戦しました。
まずは、第一段階として、最適な潜在変数パラメータを探してみました。
今回、少し苦労したのは、
-
numpyのバージョン問題
テキストマイニングの紹介記事でも書きましたが、今回も2.X系のnumpyは、動作しませんでした。
気をつけましょう。 -
PLSのパラメータ値 n_componentsの最大値問題
これはちょっと、はまってしまいました。なぜ37でエラーになるのか最初分かりませんでした。
交差検証のことも考慮しないといけませんね。勉強になりました。
次回は、今回のパラメータ値を使ったPLSモデルを実装してみます。
※ 3/4 大幅に修正し申し訳ございませんでした。モデルを作成して精度を確認した際にかなり精度が悪く
見直しました。
参考
SPSS Modeler ノードリファレンス目次
SPSS Modeler 逆引きストリーム集
SPSS funさん記事集
SPSS連載ブログバックナンバー
SPSSヒモトクブログ
今後は、ヒモトクブログなどは以下のTechXchangeのコミュニティに統合される予定です。
ご興味がある方は、ぜひiBM IDを登録して参加してみてください!!!お待ちしています。
IBM TechXchange Data Science Japan