Pirika logo
JAVA,HTML5と化学のサイト

Pirika トップ・ページ

Pirikaで化学
 物性化学
 高分子化学
 化学工学
 分子軌道
 情報化学

 その他の化学
 アカデミア
 MOOC講義資料
 プログラミング

ハンセン溶解度パラメータ(HSP):
 HSP基礎
 HSP応用
 ポリマー
 バイオ・化粧品
 環境
 物性推算
 分析
 化粧品の処方設計
 その他
 自分でやってみよう

雑記帳

Ad Space for you

 

Ad Space for you

 

 

 

Last Update

25-Jan-2019

ぴりかで化学情報化学: PLS(部分最小自乗法)

2018. 12. 29  非常勤講師:山本博志 講義補助資料

 

PLS(Partial Least Squares)法は計量経済学者であったHerman WoldとSvante Woldによって開発された新しいモデリングの方法である。主成分分析と似ているが、プログラム的には大きく異なっている。特に化合物のスペクトル・データのように説明変数が膨大になる時や説明変数間に相関がある(共線性)場合に効果を発揮する。最近のMaterials Informaticsなどでは、RDKitなどを用いて膨大な化合物情報が得られるようになった。例えば「化合物の水への溶解度を推算するのに、1000を超える化合物情報(FP: Finger Print 指紋)を用いて機械学習を行って予測を行う」などが行われている。そうしたモデリングを行う際に、PLS回帰法やニューラルネットワーク法が使われている。

PLS法自体は統計解析パッケージに搭載されているし、Pythonなどで簡単に試してみることができるが、如何せん情報が非常に限られていて手軽に試せるような状況では無いように思える。(情報が限られているという意味は、例えば統計解析用の書籍には主成分分析がほとんどでPLSは余り記載が無い。Webでは情報がありすぎて具体的にどう良いかが見えてこない。)

ブラウザー上でPLSを計算するプログラムを作ったので、それを元にPLS法の使い方を学んでみよう。

プログラムのタネ本は、「Excel VBAによる化学プログラミング、佐藤寿邦・佐藤洋子共著 培風館 2002年」である。ここに記載のVBのプログラムをJavaScriptに変換した。

慣れてきたら次の書籍で深く学べば良いだろう。
ケモメトリックス 化学パターン認識と多蛮量解析 宮下芳勝・佐々木慎一著 共立出版 平成6年
化学者のための多変量解析 ケモメトリクス入門  尾崎幸洋・字国明史・赤井俊雄/著 講談社サイエンティフィク 2002年
ケモインフォマティックス 予測と設計のための化学情報学
J.Gasteiger・T.Engel 編 船津公人 監訳
船津 公人 ・佐藤 寛子 ・増井 秀行 訳  平成17年

線形重回帰、主成分回帰、PLS回帰どれにも共通して言えることは、化合物の物性と分子構造との関係を解析したいということだろう。(場合によっては、プロセスと効率など分子では無いものも含まれる。)知りたいものを目的変数(y)と呼ぶ。知りたいものの性質を表すものを説明変数(x)と呼び、これらは複数あることがある。そこで、
y=function(X1・・Xn)
の関係を定量的に知りたい時に回帰計算を行う。ある目的変数を高くしたい(低くしたい)時に、説明変数をどのように動かしたら良いかを定量的に知りたい時と考えれば良いだろう。
(それ以外、毒があるない、燃える燃えない、甘い苦いなどを区別するのに定性的に区分するという時にも主成分分析、PLS解析は使われるが、重回帰法という時には定量的な解析を指すと自分は考えている。)
主成分回帰とPLS回帰は比較的似ている。どちらも、説明変数がスペクトルデータのように膨大であったり、実験データ数よりも説明変数の方が多いケースに使うのは有効だ。説明変数の次元数を下げるとか言う。もう一つ大事な点は、説明変数間の相関が高い(共線性のある)データが存在する場合、重回帰法では正しいモデリングができないが、主成分回帰とPLS回帰では正しくモデリングできるなどの違いがある。主成分回帰とPLS回帰の違いというのは、(素人である自分には説明が難しいが)主成分分析は逆行列を求めなくてはならない手間が大きい。行列が大きくなってくるとPLSの方が有利になる。逆にPLSには無く主成分分析にある特徴は何かあるのだろうか? 主成分同士が直交している(内積がゼロ)ことぐらいであろうか?(それが何が良いことなのかは説明できないが) ここしばらくはPLS法と重回帰法を押さえておけば良いと思う。

習うより慣れろで実際にやってみよう。
まずはプログラムがきちんと動作しているかの確認のためスペクトルデータの解析を行ってみよう。例題は種本に載っていたガソリン類の近赤外スペクトルとオクタン価の関係だ。下記のデータをコピーして表計算ソフトにペーストしておこう。

このデータは、1−36水準の近赤外のスペクトルデータと第2列目にはオクタン価が入ったテーブルになる。つまり、目的変数はオクタン価、説明変数はスペクトル値となる。
Spectrum of gasoline
全部で57種類のガソリンの36箇所での吸収強度を図に表すと上記のようになる。ほとんどのガソリンでスペクトルはほぼ一致するので、図中に57本の線があるようには思えない。この微妙な差からガソリンのオクタン価を定量的に当てると言うのが最初の例題だ。
先ほど取り出したスペクトルデータをPLS計算ソフトの上のテキストエリアにペーストする。

PLS calculation

NP値を指定してCalc. PLSボタンを押すとほぼ瞬間に計算結果が下のテキストエリアに表示される。この結果を全て選択してコピー、そして適当なテキストエディターに貼り付けておく。この出力結果の見方を覚えよう。コンピュータのプログラミングに長けているものが見たら直ぐに分かると思うが、この出力とスペクトル値があればオクタン価を予測するPLS式を簡単に組み立てることができる。プログラミング作法は別にどこかで説明しよう。ここではPLS解析結果のうち大事な項目について詳しく見ていこう。まず、ここではNP:コンポーネントの数を10と指定した。PLSの解析法は、簡単に言えば目的変数Yと説明変数Xの関係を直接求めずに潜在変数を介してモデリングを行う。最初のモデリングでYの残差がFだったとすると、その残差分をモデリングして徐々に残差を減らしていく。その繰り返し計算を何回行うかをNPとして指定する。
Residual Error of PLS

NP値が増えるにつれて、残差項の自乗和は上図のように下がってくる。NP=4であれば、ほぼ問題ない精度でモデリングできていると考えられる。

PLS regression of Octane #

そこで、NP=4で計算した訓練データの比較と、訓練に使わなかったデータの予測性能を検証した結果を上図に示した。スペクトルデータだけからオクタン価が予測できると言っても過言ではないことが分かるだろう。(ただし、Web版では予測機能は搭載していない。必要な情報は吐きだされているので、自分で組んでみよう)
具体的に言えば、PLS_W[i][j]とPLS_P[i][j]、PLS_Q[j]は出力に含まれる。PLS_Wは荷重行列(Weight Matrix), PLS_Pはローディング、PLS_Qは校正ベクトル(Calibration Vector)と呼ばれる。後は予測用のスペクトルデータPLS_XP[][]があれば次式でオクタン価予測値(PLS_YP[])が予測データセット数(PLS_MP)だけ求まる。

for(let k=1;k<=PLS_MP;k++){//予測データの数だけループを回る。
for(let j=1;j<=PLS_NP;j++){//ここで、NP回分だけループを回る。
PLS_TP[j]=0.0;
for(let i=1;i<=PLS_PW;i++){
PLS_TP[j]=PLS_TP[j]+PLS_XP[k][i]*PLS_W[i][j];
}
for(let i=1;i<=PLS_PW;i++){
PLS_XP[k][i]=PLS_XP[k][i]-PLS_TP[j]*PLS_P[i][j];
}
PLS_YP[k]=PLS_YP[k]+PLS_TP[j]*PLS_Q[j];//NP回分足し合わされる
}
}

この式の意味が理解できれば、自分で作ったPLSの結果を計算するWebアプリを幾つでもHPにおくことができると言うことになる。

石油会社では、高機能性ガソリンとして、ハイ・オクタン価のガソリン販売に力を入れている。ガソリン自体は様々なオクタン価の成分からなる混合物である。オクタン価とはn-ヘプタンを0、イソオクタン(2,2,4-Trimethylpentane)を100とした時のガソリンの燃焼効率を表す指標である。基本的には同じ炭素数の化合物なら枝分かれが多いほどオクタン価は高くなる。直鎖状の炭化水素化合物は水蒸気リフォーミングなどを通じてオクタン価をあげる操作が行われる。そうしてリフォーミングされたガソリンはどのくらいのオクタン価になったかどうしたらわかるだろうか?そのような時に、少量抜き出して近赤外スペクトルを測定し、目標とするオクタン価に達していたら反応を終了するといったプロセス管理にPLS解析法が使える。オクタン価自体は分子のトポロジーの問題としても興味深いので別途取り上げる。

このオクタン価の問題をPLSで解いた例は、佐々木先生の「ケモメトリクス」の著書の中にも取り上げられている。そこでは、43種類の無鉛ガソリンの近赤外スペクトル700点の吸光度をデジタル化して用いているので、PLSならではの例題になっている。逆に言うと、今回の例題の36ポイントでの吸光度程度では、PLSの良さを発揮できない。線形重回帰でも下図のように計算できてしまう。
Multiple Regression calc. Octane No
逆にどうすれば、700データポイントから、36データポイントを選び出すことができたのだろうか?
やり方は記載されていないが、おそらくGAPLS(遺伝的アルゴリズムを用いた変数選択PLS)などを用いて、精度を下げることなく、700個から36個を選択したのではないかと思う。PLS法は本来、共線性のデータがあっても影響を受けにくい解析方法である。しかし、NPが増えてくると影響を受け出すので、結果に影響を与えない列や共線性のあるデータはなるべく除いてから計算する方が良い。そうしたことを可能にするGAPLS法は、東京大学の船津先生や、明治大学の金子先生が得意とする解析方法である。

先ほどから共線性と言う言葉が何度か出てきている。参考図書を読めばちゃんとした定義も含め詳しく説明されているが、自分にはいまいちピンとこない。具体例で計算してしてみよう。次の例はパソコンデータ分析入門 回帰分析編(技術評論社、小幡卓+堀合啓ー著)にあった例題だ。ある石鹸工場で一月にどれだけスチームを使うか、その時の製造条件などと共に2年分が記録されている。石鹸というのは油脂を加水分解して脂肪酸とグリセリンに分離する。油脂は温度が低くなると固化してしまうのでスチームで温めなくてはならない。気温やら風速など色々なデータが収録されている。下のデータをコピーして表計算ソフトにペーストしておこう。

早速PLSを計算してみよう。

NP vs Sum of Error
NP=4あたりで精度は一定になるので、NP=4で打ち切ることにする。その時のモデルは以下のような精度になる。

Steam Usage

オクタン価の場合と比べてあまり収束性は良くない。
この問題に対して、与えられた列全てを使って重回帰法で計算をやってみよう。
MRCalc. Steam Usage
ほぼ似たような精度で、スチームの使用量が計算できる。
それでは、No1の他の条件は同じのまま、風速だけ0.2刻みで増やした時の、PLSと重回帰計算の予測値を比べてみよう。

PLS_MR Steam usage

重回帰計算(MR)では風速が早くなるとスチームの使用量は減るのに、PLS回帰ではわずかであるがスチームの使用量は増えると予測する。どうしてこのような違いが起きるのだろうか? これが共線性の問題である。

Steam usage vs Temp
もともと、気温が下がるとスチームの使用量は増えることは上図を見なくても納得できるだろう。


WInd Temp
しかし、当たり前の事ながら、冬場は気温が低く風速も早い。このように説明変数間に相関があるような場合を「共線性がある」と言う。
重回帰式の係数は次のようになる。
スチームの使用量=・・・・・・-0.089167839816838*風速
-0.0825515068788824*気温+・・・・+4.38023018565945
風速の係数は符号がマイナスなので風速が早くなるとスチームの使用量はより減ってしまう。通常であれば、風速が早い時は(冬場で)気温が低い。気温の重回帰係数もマイナスなので、気温が低いと言うことはマイナス分が小さくなるのでスチームの使用量は増える。つまり、風速が早いということと気温が低いということは逆相関であって共線性がある。そこで重回帰係数は不安定になるので正しくスチーム使用量を予測できなくなる。脂肪酸量とグセリン量にも正の相関が存在する。このような問題に対してはPLS法を使った方が解の安定性は高まるとされている。

ただし、この程度の説明変数量であるなら、変数選択重回帰法を取れば良い。
3Param Steam Predic.
スチームの使用量に対して8説明変数が上げられていたが、実際には3変数を選択しただけでも上図のように相関が取れてしまう。3変数まで減らしてしまえば、風速と気温のような共線性のある項目は選ばれることはない。予測式を作る際の多くの間違いは精度を上げすぎるところにあると思う。NP値を上げすぎる、説明変数を使いすぎる。見かけは精度が高くなったように見えるかもしれないが、Trainデータに対する記述性が高くなることと、Testデータに対する予測性能が高いこととは別問題である。PLS法に頼るだけでなく、説明変数の意味をきちんと理解して取捨選択する化学者の能力が大事なのは言うまでもない。

RDKitとPLS法を用いた成果のうち最も著名なのは化合物の水への溶解度の例であろう。RDKitをダウンロードするとTrainとTestのsdf・データが含まれているのですぐに試してみることができる。Pythonから扱うのであれば、sdfファイルを直接変換すれば良いだろう。Pirikaのソフトウエアーを利用するなら、solubility.test.sdfは気をつける必要がある。このテキストファイルは改行コードが\r\nになっていることと、ファイルの途中でSmiles, SMILESの大文字小文字が混在するので、少し修正しないとSmilesのテーブルは得られない。しかし、一旦Smilesのテーブルさえできてしまえば、それからRDKitを使って識別子を吐き出させ、PLSで予測式を作成し、以下のような図を作成するのには、5分とかからない。 2019.1.9:分子結合インデックスの計算に誤りがありました。本来、Chiパラメータは水素原子は計算に含めないはずです。Smilesの構造から3次元分子を作る際に水素を付加させていました。RDKitでは水素の有る無しで計算結果が変わってしまいます。

似たような検討は大学などでも行なっています

Solubility to water calc. by PLS

大学での計算例は、RDKitのFP:Finger Printなどを使っているようであるが、結果はほとんど違いがないように思える。(RDKitの吐き出す識別子はPro版を使っているので117種、PLSの計算もPro版を使っているので予測値も簡単に計算できるが、必要な情報は吐き出されているので、誰でも同じことは検証できるだろう。企業の研究者など手間を節約したい場合は、Pro版の利用を考えて頂きたい。)水への溶解度(g/100cc水)のlogをとったものなので、予測値が2違うという事は、100倍溶解度が異なるということで、これを見て十分な精度と言えるかどうかは議論の余地が残るだろう。さらにデータ数が増えた時にどうなるか?こちらで解説を行なっている

こうした例題を通じてさらに深くPLSを使ってみたいのなら、先に上げた参考文献をよく読むことをお勧めする。

PLS法は新しい解析方法なので、どんな問題に適用しても論文化できると間違った風潮があるような気がする。スペクトル解析など説明変数が非常に多い系を除けば、大抵の問題は、やはり30年前の古い技術、変数選択重回帰法でも解けてしまう。変数選択の時に共線性のデータを除外するプログラムなど非常に簡単に書けるだろう。

SMILESの構造式と水への溶解度のペアがたくさんあったとする。AIが勝手にRDKitで計算し、分子軌道計算し、自動的に説明変数を集める。ネット上にある全ての水への溶解度のデータに適合する予測式をPLSや変数選択重回帰法を使ってAIが作成してしまう。そんな時代は来るのだろうか?
コンピュータの中で自己対戦を繰り返し、目的とする数式の性能アップを計るのはコンピュータの一番得意とする分野だ。特にPLSを使ったものは自動化しやすいので人間の研究者もうかうかはしていられない。MAGICIAN養成講座の大事なポイントだ。