化学研究で役に立つデータ解析入門：回帰分析の活用を広げる編

[スポンサーリンク]

一般的な話題

化学研究で役に立つデータ解析入門：回帰分析の活用を広げる編

[スポンサーリンク]

前回の化学研究で役に立つデータ解析入門：回帰分析の応用編では、Rを使ってエクセルにはできない回帰分析を行いましたが、データセットが大きくなってくると、どの変数を取り込むべきか悩むことがあります。そこで今回は変数を選ぶ際に有用な機能を紹介します。

回帰分析の活用を広げる方法

具体的にステップワイズ、多重共線性、統計量について本記事では取り上げます。

ステップワイズ：有機合成の実験において反応時間や触媒量、反応温度などが収率に影響を及ぼすと一般的に知られています。しかし、基質や溶媒、触媒の純度などの物質に関するパラメーター、実験室の温度・湿度などの環境に関するパラメーター、スターラーバーやフラスコの形・清浄性といった実験器具のパラメーター、実験者や分析装置といった作業に関するパラメーターなど、実は収率に影響を与える要因はたくさんあり、通常はそれらをなるべく合わせて実験を行いますが、全てのパラメーター揃えることは現実的ではありません。何度か同じ実験を繰り返しても収率が変化する場合には、何か隠れた条件に対応して収率が変化していると言えます。収率を変化させているパラメーターを探す時にも回帰分析を活用できますが、たくさんの変数を取り込むと、R2は高くなるもののP値も高い変数ばかりでどの変数の影響が大きいのか判断が難しくなります。そのため有効な変数を選択してモデルに組み込む必要があり、ステップワイズという機能を使うと有効だと思われる変数を自動的にピックアップすることができます。

多重共線性：コンビニの夏のアイスクリームの売り上げについて調べたところ最高気温と天気に相関があることがわかりました。しかしながら気温と天気には相関がある（夏は天気が良いほど温度は上がる）ため、正しい解析とは言えません。ここまではパラメーター同士に関係がないことを前提にして解析を行ってきましたが、たくさんの変数をモデルに組み込むとお互いに関係がある複数の変数を含んでいる可能性が高くなります。そこで変数同士に関係性があるかないか（多重共線性）を調べることで関係が重複している変数を排除することできます。

回帰診断図：回帰分析を活用した結果をまとめて上司や教授に報告したら、「私は、（直感的に）直線の傾きがもっと緩やかだと思う」や「過去の経験ではそんな相関はない」と言われて否定されるかもしれません。そんなときには回帰の診断が必要であり、これにより結果の根拠を知ることができます。

ステップワイズ

実際にRを使ってステップワイズを使った回帰分析を行います。追加するコードは一行だけで、ステップワイズを行いたい回帰分析をStep()に入れるだけです。

回帰ファイル名<- lm( 目的変数 ~ 説明変数, data=データの名前)

結果ファイル名<- step(回帰ファイル名)

例えば、ニュージャージー州の有料高速道路での車の速度超過を調べたデータセットEGViolatorsを使うと下記のようなコードになります。

ステップワイズを行うコード

fit1とres1は任意の名前であり、Overlimitという目的変数に対してExitやDirection, Speed, Licenseを説明変数を代入して相関があるか調べました。データセットにはRace=人種も記録されてますが、数字で表記されていてどの数字がどの人種を示しているか不明だったので省きました。分析を実行すると下記のような結果が得られます。

ステップワイズの過程

RではAIC(Akaike’s Information Criterion)という方法で説明変数の選択が行われ、、ざっくり言えば適度に有意義な変数を適度な数だけピックアップされます。各ステップごとに次にどの変数を抜くとAICの値がどうなるかを示していて、一番下では、<none>＝何も変数を抜かないが最小値を示したので、ExitとSpeedで最適解が得られたと判断されました。

ステップワイズ後の線形回帰

結果、P値が極めて低い変数が見つかり、Exit番号が高く、スピードが速い方が速度超過も大きいことが分かりました。（高速道路全体で制限速度が異なる前提のため、スピードが遅くても制限速度が低いと速度超過が大きいこともありうる考えられます。）

このようにステップワイズを使うとたくさんの変数の中から有効だと思われる変数を絞ることができます。もちろん分類変数に対しても有効なので、合成実験において再現性の問題が出た場合にはなるべくたくさんの項目（天気や湿度、昼食のカロリー、~~先輩の機嫌~~）を入れてこの分析を行うと、原因がわかるかもしれません。

多重共線性

多重共線性は、Variance Inflation Factor (VIF)を計算することで評価できます。Rでは、carパッケージのvif()を使うと簡単に計算できます。

library(car)

回帰ファイル名<- lm( 目的変数 ~ 説明変数, data=データの名前)

vif(回帰ファイル名)

VIFを算出するコード

carパッケージを使うために冒頭には常にlibrary(car)を入れ、vif()に調べたい回帰分析名を入れるだけです。エクセルでのデータ解析のときに使用したねじの数を数える機械についてのデータセットを使い、あえてすべての変数を入れてvifを計算します。

VIFの結果

結果、TIMEとTOTALが高い結果になりました。VIFが高いほど多重共線性が疑われ、TIMEとTOTALには関連があると言えます。加えたねじが多いほど機械が全てのねじを数えるのに時間がかかるので、実験事実からもこのVIFが高い理由は説明でき、20個のねじを数えるのにかかる時間を予測するためには両方の変数を使うべきではないと言えます。

このねじの場合は単純で、回帰分析を行う前に変数の関係性を把握できますが、変数が多くなってくると関係性を掴むことは難しくなります。またステップワイズではこの多重共線性は考慮されないため、VIFが高い変数が選ばれてしまう可能性もあります。vif()は分類変数でもステップワイズの結果に対しても使えるので、常にチェックしておくべきだと思います。一般的にはVIFが10を超えると多重共線性があるとの解釈になります。

回帰診断図

いろいろな診断がありますが、plot()を使うと，「残差と予測値の散布図」，「残差のnormal Q-Q プロット」，「標準化した残差の絶対値の平方根と予測値の散布図」，「Cookの距離」が表示されます。

使い方は簡単で、plot()に調べたい回帰分析名を入れるだけです。Rstudioでは、ConsoleにReturnを押してくださいと表示され、Enterを押すとグラフ表示エリアに4つの図が順々に表示されます。左矢印で前の図に戻ることができ、Exportでお好みの形式で保存もできます。ここでは、ねじの数を数える機械についてのデータセットを使い、T20BOLTに対してRUNとTOTAL、SPEED1を説明変数としました。

ねじのデータセットにおけるResiduals vs Fitted

最初のグラフ、残差と予測値の散布図は、得られた回帰曲線（この場合はT20BOLT=0.43XRUN+0.71XTOTAL+10.65XSPEED1-31.69）と実測値のずれを表すグラフです。X軸が目的変数の値（この場合はT20BOLT）で、Y軸が残差＝予測値とのずれ（秒）を示し、各プロットが真ん中に近いほど、実測と予測が近いことになります。予測よりかけ離れた値は、データポイントの列番号が表示され（19, 20, 24）表にてデータを確認することができます。T20BOLTが30から55付近は、残差が大きく何か理由があるかもしれません。

ねじのデータセットにおけるnormal Q-Q plot

２番目に表示されるグラフ、残差のnormal Q-Q プロットは、得られたデータと理論分布を比較しその類似度を調べるためのグラフです。X軸に正規分布に従う場合の期待値を横軸にとり、Y軸の値は規準化した残差をとります。正規分布に従う場合には、きれいにプロットが点線に乗ることになります。Residuals vs Fittedと異なり予測値と実測値の振れが正規分布かどうかを調べる図と言えます。予測よりかけ離れた値はデータの列番号が表示され、このデータではResiduals vs Fittedと同様に19, 20, 24番が振れが大きいところで正規分布からずれていることを示しています。

ねじのデータセットにおけるscale_location

３番目に表示されるグラフ、標準化した残差の絶対値の平方根と予測値の散布図もResiduals vs Fittedと同様に得られた回帰曲線と実測値のずれを表すグラフです。ただし、Y軸の値は規準化した残差の絶対値の平方根となっていて、残差の変動状況を考察するために使用されます。

ねじのデータセットにおけるResiduals vs Leverage

最後のグラフ、Cookの距離は、X軸にある実測値からデータセット内のすべての実測値の平均までの距離（てこ比）を示し、Y軸の値は規準化した残差をとります。同時にこの図の中に、クックの距離が0.5と1の範囲を赤の点線で示されます。Cookの距離は、個々のデータが回帰式の推定に及ぼす影響を表した距離のことで、値が大きいほど回帰式の推定に大きく影響していると言えます。このデータでは、クックの距離が大きいプロットはなく、線は図の外なので表示されていません。

プロット別のクックの距離自体を表示させるには下記のコードを使用します。

plot(cooks.distance(回帰ファイル))

X軸がデータの列番号でY軸がクックの距離を示し、このデータセットでは最大でも0.1ほどで際立って高いデータはなく、特定のデータが回帰式を大きく傾けていることはないとわかります。

ねじのデータセットにおけるクックの距離

このように、最初に表示されるグラフほど直感的に理解できるデータで、後に表示されるグラフほど回帰式への影響を直接的に示すようなグラフとなっています。集団から大きく外れたデータポイントがなければ、まんべんなく影響を受けたモデルだと主張することができ、逆にクックの距離が異常に大きいプロットがResiduals vs Leverageで見つかれば、振れの符号と大小を他のグラフから調べて、そのデータポイントを除くか補正するかを検討できます。本記事ではそれぞれの例で必要なデータしか出力していませんが、いろいろな結果が必要な場合には出力するコードを羅列すれば、一度にいろいろな種類のデータを取得することができます。

データポイントと変数が多くなるほど、解析は複雑になりいろいろなパラメーターを出し入れして最適なモデルを作る必要があります。そんなときに、ステップワイズで機械的に選択し多重共線性を確認したり、回帰診断図でどんなデータポイントが影響を与えているか確認することがより当てはまりの良いモデルを作るのには有用になると思います。また各値の数学的な説明は省きましたが、導出方法を知ることでよりその値を意味を理解することができます。