[スポンサーリンク]

一般的な話題

化学研究で役に立つデータ解析入門:回帰分析の応用編

[スポンサーリンク]

化学研究で役に立つデータ解析入門:エクセルでも立派な解析ができるぞ編ではエクセルでも立派な回帰分析ができることを紹介しました。しかし、回帰分析は奥が深く、エクセルの機能では、できないこと、できても手間がかかることがたくさんあります。今回は、フリーソフトウェアの統計解析向けのプログラミング言語Rを使ってエクセルよりも簡単により深く回帰分析を使いこなす方法を紹介します。

エクセルの回帰分析ではできないこと

回帰分析は奥が深く、エクセルではできない使い方や計算されない統計量がたくさんあります。化学実験の結果を解析するにあたっては下記の機能が特に有用であり、今回はこれらに焦点を当てて解析を行っていきます。

分類変数の取り扱い:実験の条件は、温度や反応時間といった連続した数字(連続変数)だけでなく、基質の官能基や触媒の種類といった数で表すことができない条件(カテゴリー変数)もあり、これらの傾向を掴むことも重要な解析となります。しかしエクセルでカテゴリー変数を取り扱うには、各種類(例えば官能基)ごとに1か0で表現(ダミー変数)して解析する必要があり、解析前のデータの加工が必要になります。

信頼区間:信頼区間とは母集団の真の値が含まれることがかなり確信できる数値範囲のことで、この”かなり確信できる”は、ばらつきの状況によって変わり、95%(95%信頼区間)が一般的です。具体的を95%信頼区間を表すと、ある条件で100回実験を行ったら95回の結果が含まれる結果の範囲であり、予測の誤差を直感的に知るのに便利な統計量です。エクセルでも各係数における信頼区間は算出されていますが、結果を予測する際に重要なのは任意の条件における信頼区間であり、それは画像検索するとたくさんヒットするラッパ型の曲線を算出することで、あるxにおいてyがどれくらいの幅があるのかを知ることができます。エクセルの回帰分析で算出される統計量でもグラフを作ることは可能ですが、手作業でグラフを作る必要があります。

交互作用:例えばある合成において反応温度を上げれば収率は単調増加するが、同時に触媒量を増やすことで収率の増加量が加速する場合、反応温度と触媒量に正の相互作用があると言えます。逆に反応温度を上げれば収率は単調増加するが、同時に基質の量を増やすと収率の増加量が抑えられる場合。反応温度と基質量に負の相互作用があると言えます。エクセルで相互作用を調べるには単純に、条件同士を掛け合わせた新しい条件を追加して回帰分析を行えば良いのですが、相互作用ごとに行を追加すること必要で、条件が多い場合には不向きです。

Rとは

統計解析に特化したソフトウェアはいくつかありますが、ここでは、フリーのプログラミング言語Rを使っていきます。プログラミング言語といわれると、テキストエディタを使って真っ白なところにガリガリとコードを書いて何度も実行してやっとほしい結果が得られるイメージを自分は持っていますが、このRは非常に簡略的にコードを書くことが可能で、前回行ったボルトの数を数える回帰分析解析もたった2行で同様の結果を得ることができます。

前回行ったエクセルで行った解析を行うためのコード

Rstudioと呼ばれる統合開発環境もフリーで提供されていて、コードの記述支援やデータの確認も容易にできるようになっています。

Rstudioのインターフェイス

RとRstudioインストールの仕方は、ネット上にたくさん紹介されていますが、一例として下記のサイトを挙げます。

Rによる回帰分析の基礎

有用な解析に移る前にRでの回帰分析の基礎的な方法を紹介しますが、この記事では、どのようにコードを記述すれば解析ができるかについて注力し、記号の意味などは取り扱いません。コードのルールなどを知りたい場合には他のページを参照しください。また、コードは同じ機能でもいろいろな書き方があり、他のサイトでは異なるコードを使っているかもしれません。

まず、データを読み込むために右上のImport Datasetをクリックして目的のファイルを選択します。csvやtxtの場合は、From Text(base)か(readr)で、エクセルファイルの場合にはFrom Excelから読み込むことができます。エクセルで他形式のデータを読み込むときと同様に区切りの識別によって配列を確定させることができます。読み込みも x <- read.table(“data01.txt”)などを使えばスクリプトでデータを読み込むことも可能ですが、Import Datasetを使うとデータが正しく読み込まれているか確認できるので確実性があります。読み込まれたデータセットには名前(拡張子抜きのファイル名)が付けられて、ダブルクリックすると内容を確認することができます。

各ウィンドウの機能(表示されているタブの機能であり、タブを変えると異なる機能となります。)

回帰分析は、lmというコードを使います。最低限の記述で表すと、

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

で回帰分析を行うことができます。

  • 結果ファイル:結果を格納するデータの名前、何でもよい
  • 目的変数:結果の項目(エクセルでいうY範囲)で読み込む表のカラムの名前を入れる
  • 説明変数:変数の項目(エクセルでいうX範囲)で目的変数同様に読み込む表のカラムの名前を入れる。複数の変数を取り込む場合には+で変数を足していく
  • データの名前:調べるデータの名前で、読み込んだデータの名前を入力する

基本の解析はこれだけです。この一行を実行する(Sourceをクリック)と、解析が行われて結果が入ったデータが作られます。表を確認するときと同様にデータをダブルクリックしても解析結果を確認できますが、もう1行summary(データの名前)と記述するとConsoleに統計量が表示されます。

実際のコード

上記のコードを実行した時の結果、エクセルの結果と全く同じことがわかる

分類変数がある回帰分析

では、エクセルで処理できない解析方法をやってみます。分類変数を取り扱う場合も方法は同じです。ボルトを数えるデータセットには分類変数がなかったので、新たなデータセットEGViolatorsを使います。これは、ニュージャージー州の有料高速道路での車の速度超過を調べたもので、黒人の方がスピードを超過しているという疑惑を解明するために超過速度を人種や運転免許の発行州、測定地点、走行方向などとともに調べた結果です。黒人と他の人種で超過速度に差がないことが証明されていますが、今回は目的を変えて、運転免許の発行州別に速度超過の傾向に差があるのかを調べます。

実際に回帰分析を行う方法ですが、上記と全く同じで、運転免許の発行州であるLicenseを代入するだけです。

OverlimitとLicenseの関係を分析するコード

実行すると運転免許ごとに値が算出されます。これはAL(アラバマ州)の運転免許を基準に他の州の運転免許では超過速度がどう変化したを示しています。式で表すなら超過する速度=5.67(InterceptのEstimate)+免許のEstimateとなり、例えば、表の真ん中付近のMN(ミネソタ州)の運転免許保有者は5.67+7.83=13.5 マイル/hほど統計的に有意差がある中で速度超過していることになります。発行州によっては、P値が高く有意差がないと言える場合もありますが、ミネソタ州のように有意差がある場合もあります。ただし、Adjusted R-squared(補正R2): 0.005996 が1に程遠いので、全体としては速度超過と運転免許の発行州は相関が小さいと言えます。

回帰分析の結果、P値が低さによって右端に記号が付与される。

このようにRでは分類変数を組み込んでも問題なく解析されます。結果は値の場合とは異なり、ある基準の分類(上記ではAL)と比べた係数の大小で表されます。

信頼区間の図示

そもそもエクセルでは、回帰分析で得られた結果を実測値をプロットする機能はなく、手作業で係数を拾ってグラフを作る必要があります。しかしRを使えば5行追加するだけでグラフまで作ることができます。ここでは新たなデータセット、1970年代後半におけるボストンの住宅価格をデータセットとして使います。これは様々な条件と町の住宅価格の間に相関があるか調べたものです。

まず散布図を作成するには、

plot(Y値~X値, data=データ名)

で、lmと同様にXとYに代入したい表のカラムを入力し、データ名に取り込んだデータセットの名前を記述するだけです。例えば、AGE:1940年より前に建てられた持ち家の割合=「古い家の割合」をXとしてMEDV:「住宅価格」(1000ドル単位)の中央値でプロットすると図が右下に表示されます。

X軸にAGE、Y軸にMEDVをプロットした散布図を作成するためのコード

グラフを作成した結果

とてもシンプルでわかりやすいですが、タイトルを入れたりプロットの色を変える場合には、コードを追加してグラフをデコレーションすることができます。

では、回帰分析を行って曲線を追加します。追加するコードは、

abline(回帰分析結果)

で、先ほどのボストンの住宅価格の例を使うと、回帰分析後にablineを入れるだけです。plotによって作られた散布図に追記する形をとるので、このコードより前に散布図が作られていないとエラーになります。

回帰分析を行ってその曲線をグラフに上書きするコード、col=’red’で曲線の色を赤色にしている。

直線を追加したグラフ

では最後に信頼区間を表示させます。

配列名<- data.frame(X軸に使うカラム名 = seq(始点, 終点, 感覚))
信頼区間名 <- predict(回帰分析結果, newdata = 配列名, interval = ‘confidence’, level = 信頼区間の%)
lines(配列名$X軸に使うカラム名, 信頼区間名[, 1])
lines(配列名$X軸に使うカラム名, 信頼区間名[, 2])
lines(配列名$X軸に使うカラム名, 信頼区間名[, 3])

信頼区間を表示させるには、範囲の設定が必要でそれを一行目で行います。次の行で信頼区間を計算し、それをもとに3~5行目で曲線を描きます。[,1]が回帰曲線そのもので[,2]が下限、[,3]が上限に相当します。同様にボストンの住宅価格の例を使うと下記のようになります。

回帰分析後、信頼区間を計算しグラフにプロットするコード

95%信頼区間を追加したグラフ、任意のAGEにおいてMEDVの総結果の95%が緑の間に入ることを示している。

信頼区間の違いを確認するためにZN: 25,000平方フィートを超える区画に分類される住宅地の割合=「広い家の割合」を使ってプロットしました。

ZNとMEDVの信頼区間を示したグラフ

確かに広い家の割合が多いほど町の住宅価格は上昇しますが、ばらつきも大きくなってZNが80%以上では、95%信頼区間が5000ドル以上に広がっているのがわかります。商品の有効性を示す近似曲線は、日常生活でも多々見かけますが、効果があるように見えても、ばらつきが大きく実際は有効性が疑わしい場合もあると思います。そのため統計的な信頼区間の議論はばらつきが大きい結果ほど重要であると言えます。

もっと簡単に、モダンなグラフを書きたい場合には、ggplot2という追加パッケージが便利です。Rではパッケージを追加することで機能を追加することかでき、具体的には、右下のウィンドウの「Packages」のタブで追加したパッケージ名にチェックを入れるだけ使えるようになります。先ほどのAGEとMEDVの値を使って信頼区間の入ったグラフを作るには、

ggplot2を使って信頼区間を図示するためのコード

の2行で下のようなグラフを作ることができます。

ggplot2を使った作成したグラフ

もちろん、オプションのコードを増やすことでグラフの体裁を変えることもできます。

相互作用がある回帰分析

最後に相互作用がある回帰分析の方法について紹介します。相互作用を調べるときにはX値の書き方を変えるだけです。速度超過のデータセットを使って、超過速度と測定した出口(Exit)と運転免許の発行州(Licence)について関係があるかどうか回帰分析を行う場合には、result <- lm(Overlimit~Exit+License, data=Violators)と記述すると、下記のような結果が得られます。相互作用を考えなければ、超過速度と測定した出口の相関が大きいものの、特定の運転免許の発行州でも有意差がみられています。

ExitとLicenseの回帰分析の結果、補正R2は0.01034である。

そこで、相互作用を考慮した回帰分析を行います。記述をresult <- lm(Overlimit~Exit*License, data=Violatorsとするだけで、Overlimit=a×Exit+b×License+c×Exit×License+dの式の近似を行います。この結果は下記のようになり、運転免許の発行州の有意差が見られるようになった一方、Exitの有意差は見られなくなりました。相互作用の効果を見ると多くの州において負を示しています。この解釈を考えると、ExitのEstimateは正なので、番号が大きいExitほど速度超過が大きいが、多くの特定の発行州の免許保有者はその比例を強めるか弱めることになります。例えば、CO(コロラド州)の免許保有者は、平均よりも超過速度は高いが、番号が大きいExitほど、超過速度はより低くなる傾向があるということになります。全体のモデルの当てはまりを示す補正R2の相互作用を考慮することで若干向上しました。

相互作用を考慮したExitとLicenseの回帰分析の主成分の結果

相互作用を考慮したExitとLicenseの回帰分析の相互作用の結果(一部のみ掲載)

相互作用を考慮したExitとLicenseの回帰分析の結果、補正R2は0.01114である。

三つ以上の変数(例えばAとBとC)を*で接続して記述すると、A,B,C,AB,AC,AC,ABCとすべての相互作用が計算されるので特定の相互作用を除きたい場合には、考慮したい組み合わせを+で合わせます。

相互作用を考慮して当てはまりが良くなるかどうかは、データセット次第であり、当てはまりが悪くなるか条件の有意差が低くなることもあります。また相互作用を算出にはそれなりのデータ数も必要で、速度超過のデータセットでも発行州によっては相互作用の値がNAになっているのは、データ数が少なすぎて算出できなかったためです。そのため相互作用まで考慮する意味があるかどうかは、相互作用ある場合とない場合の回帰分析を比較して判断する必要があります。このように相互作用がある系の回帰分析もRでは簡単に分析することができ、単体の効果は当たり前だが、複数の条件の組み合わせによって予想外の結果が得られることはよくあり事象で、それを統計的に探して発見できることが、この回帰分析の面白さだと思います。

今回は、Rを使った回帰分析を紹介してきました。Rは世界的に認められたソフトウェアであり、FDAに薬のデータを提出する際にも使われています。また、いろいろなパッケージが開発されていて、Rで最新の解析方法を取り扱うこともできます。次回は、より深く回帰分析を使いこなす方法を紹介します。

関連書籍

関連リンク

Zeolinite

投稿者の記事一覧

ただの会社員です。某企業で化学製品の商品開発に携わっています。社内でのデータサイエンスの普及とDX促進が個人的な野望です。

関連記事

  1. マテリアルズ・インフォマティクスにおけるデータの前処理-データ整…
  2. 計算化学者は見下されているのか? Part 1
  3. 提唱から60年。温和な条件下で反芳香族イソフロリンの合成に成功
  4. ポンコツ博士の海外奮闘録⑫ 〜博士,今と昔を考える〜
  5. 立体規則性および配列を制御した新しい高分子合成法
  6. Wiley社の本が10%割引キャンペーン中~Amazon~
  7. 今年は共有結合100周年ールイスの構造式の物語
  8. 製薬業界の現状

コメント、感想はこちらへ

注目情報

ピックアップ記事

  1. 第94回日本化学会付設展示会ケムステキャンペーン!Part II
  2. ノーベル化学賞2011候補者一覧まとめ
  3. 複数のイオン電流を示す人工イオンチャネルの開発
  4. マテリアルズ・インフォマティクスの普及に取り組む事業開発ポジションとは?〜最前線メンバー3人のキャリア選択と今について語る〜
  5. 単分子レベルでの金属―分子接合界面構造の解明
  6. 2007年文化勲章・文化功労者決定
  7. 低分子医薬に代わり抗体医薬がトップに?
  8. ライセルト反応 Reissert Reaction
  9. 高分子鎖の「伸長」と「結晶化」が進行する度合いを蛍光イメージングで同時並列的に追跡する手法を開発
  10. アルデヒドのC-Hクロスカップリングによるケトン合成

関連商品

ケムステYoutube

ケムステSlack

月別アーカイブ

2020年9月
 123456
78910111213
14151617181920
21222324252627
282930  

注目情報

最新記事

Christoper Uyeda教授の講演を聴講してみた

bergです。この度は2024年5月13日(月)に東京大学 本郷キャンパス(薬学部)にて開催された「…

有機合成化学協会誌2024年5月号:「分子設計・編集・合成科学のイノベーション」特集号

有機合成化学協会が発行する有機合成化学協会誌、2024年5月号がオンライン公開されています。…

電子のスピンに基づく新しい「異性体」を提唱―スピン状態を色で見分けられる分子を創製―

第614回のスポットライトリサーチは、京都大学大学院工学研究科(松田研究室)の清水大貴 助教にお願い…

Wei-Yu Lin教授の講演を聴講してみた

bergです。この度は2024年5月13日(月)に東京大学 本郷キャンパス(薬学部)にて開催されたW…

【26卒】太陽HD研究開発 1day仕事体験

太陽HDでの研究開発職を体感してみませんか?私たちの研究活動についてより近くで体験していただく場…

カルベン転移反応 ~フラスコ内での反応を生体内へ~

有機化学を履修したことのある方は、ほとんど全員と言っても過言でもないほどカルベンについて教科書で習っ…

ナノ学会 第22回大会 付設展示会ケムステキャンペーン

ナノ学会の第22回大会が東北大学青葉山新キャンパスにて開催されます。協賛団体であるACS(ア…

【酵素模倣】酸素ガスを用いた MOF 内での高スピン鉄(IV)オキソの発生

Long らは酸素分子を酸化剤に用いて酵素を模倣した反応活性種を金属-有機構造体中に発生させ、C-H…

【書評】奇跡の薬 16 の物語 ペニシリンからリアップ、バイアグラ、新型コロナワクチンまで

ペニシリンはたまたま混入したアオカビから発見された──だけではない.薬の…

MEDCHEM NEWS 33-2 号「2022年度医薬化学部会賞」

日本薬学会 医薬化学部会の部会誌 MEDCHEM NEWS より、新たにオープン…

実験器具・用品を試してみたシリーズ

スポットライトリサーチムービー

PAGE TOP